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Summary 



Nucleons, i.e., protons and neutrons, are composed of quarks and gluons, whose interactions 
are described by the theory of quantum chromodynamics (QCD), part of the standard model 
of particle physics. This work applies lattice QCD to compute quark momentum distributions 
in the nucleon. The calculations make use of lattice data generated on supercomputers that 
has already been successfully employed in lattice studies of spatial quark distributions (" nu- 
cleon tomography"). In order to be able to analyze transverse momentum dependent parton 
distribution functions, this thesis explores a novel approach based on non-local operators. 
One interesting observation is that the transverse momentum dependent density of polarized 
quarks in a polarized nucleon is visibly deformed. A more elaborate operator geometry is 
required to enable a quantitative comparison to high energy scattering experiments. First 
steps in this direction are encouraging. 



Zusammenfassung 

Nukleonen, also Protonen und Neutronen, bestehen aus Quarks und Gluonen, deren Wechsel- 
wirkung durch die Quantenchromodynamik (QCD) innerhalb des St andardmo dells der Teil- 
chenphysik beschrieben wird. Diese Arbeit nutzt Gitter-QCD zur Berechnung von Quark- 
Impulsverteilungen im Nukleon. Dabei wird auf Gitterdaten zuriickgegriffen, die auf Hochleis- 
tungsrechnern erstellt und bereits erfolgreich fur die Analyse der raumlichen Quarkverteilung 
(" Nukleontomographie" ) eingesetzt wurden. Fur die Untersuchung von Transversalimpuls- 
abhangigen Partonverteilungsfunktionen stellt diese Arbeit ein neuartiges Verfahren, basie- 
rend auf nicht-lokalen Operatoren, vor. Die Ergebnisse zeigen u.a. eine sichtbare Verformung 
der Transversalimpuls-abhangigen Dichte polarisierter Quarks im polarisierten Nukleon. Ein 
quantitativer Vergleich mit Streuexperimenten erfordert kompliziertere Operatorgeometrien. 
Erste Schritte in diese Richtung sind vielversprechend. 
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Chapter 1 

Introduction 



What does a proton look like inside? Today, it is well established that nucleons, i.e., protons 
and neutrons, are composed of quarks and gluons, whose interactions are described by the 
theory of quantum chromodynamics (QCD), part of the standard model of particle physics. 
Quarks carry a charge called "color" , and are kept confined in color-neutral bound states by 
the strong interaction mediated by the gluons. Consequently, we can never observe isolated 
quarks in a detector. The proton has a radius of the order of lfm = 1 x lCP 15 m, inside of 
which the quarks are confined. According to Heisenberg's uncertainty principle, this means 
that the quarks cannot be permanently at rest in the nucleon. There is an intrinsic motion 
of quarks inside the nucleon, or, more precisely, the quark momenta follow distributions of 
non-zero width. In this work, we investigate this momentum distribution with theoretical 
means, using lattice QCD. 

Experimental information about the internal structure of the nu- 
cleon is obtained in high energy scattering experiments, in which we 
can probe the nucleon at large velocities close to the speed of light. 
Due to the Lorentz contraction, such a fast nucleon appears flat 
like a disk, as illustrated in Fig. 4.1. During the scattering exper- 
iment, a hard interaction takes place with one of the quarks inside 
the nucleon. At the instant of the interaction, the quark carries 
a momentum fraction x of the nucleon momentum P, and has an 
intrinsic transverse momentum k± perpendicular to the direction 
of flight of the nucleon. The spatial location of the quark in the 
transverse plane is given by a vector b±, the so-called impact pa- 
rameter. The distribution of quarks with respect to the longitudinal 
momentum fraction x is experimentally most easily accessible, e.g. 
inelastic scattering (DIS) experiments. It is parametrized in terms of conventional parton 
distribution functions (PDFs), such as fi(x). In general, PDFs tell us how likely it is to find 
a quark that carries a given momentum fraction x of the nucleon. The concept of PDFs can 
be extended to transverse momentum dependent parton distribution functions (TMD PDFs) 
such as fi(x, fejjj which describe quark distributions with respect to both longitudinal mo- 




Figure 1.1: Illustration 
of a nucleon with large 
momentum P. 

in fully inclusive deeply 
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mentum x and transverse momentum k±. Effects of the intrinsic transverse quark momentum 
are observable, e.g., in the angular distribution of the measured final state particle in semi- 
inclusive deeply inelastic scattering (SIDIS) experiments. Another extension of PDFs are the 
generalized parton distribution functions (GPDs). Those enable us, amongst other things, to 
visualize the nucleon in terms of &j_ -dependent quark densities [BurOO, RP02, Bur03, DH05] 
( "nucleon tomography" ) . 

In the study at hand we attempt to calculate transverse momentum dependent parton distri- 
butions from first principles, i.e., from the theory of quantum chromodynamics. Although the 
laws of QCD are given by a few elegant mathematical expressions, it is very difficult to cal- 
culate properties of nucleons or other hadrons, due to the strong interaction. In recent years, 
great progress has been made within lattice QCD, a method that allows us to perform quantum 
field theoretical "simulations", albeit at an extreme computational cost. In particular, it has 
become possible to calculate GPDs [G+04, H+04, S+04, R+05, G+05b, D+05, G+07, H+08], 
which yield "tomographic" pictures of the nucleon, showing an interesting deformation in the 
spin sensitive channels, see, e.g., Ref. [G + 07]. 

The goal of this work is to devise methods within lattice QCD that enable us to obtain similar 
results for TMD PDFs. The fcj_-dependent distributions provide information complementary 
to the bj_-dependence encoded in GPDs 1 . The ability to visualize the fcj_-dependence of 
quark distributions inside the nucleon has been an important incentive for our study of TMD 
PDFs on the lattice. In contrast to experiment, it is particularly easy on the lattice to study 
spin dependent effects. Analyzing polarized quarks in a polarized nucleon, we shall find an 
analogous deformation as seen, previously, with nucleon "tomography" in the bj_-plane. 

The motivation to study TMD PDFs is not solely based on our interest in the nucleon it- 
self. In many scattering experiments, including those exploring physics beyond the standard 
model at the LHC, it is neccesary to give an accurate description of hadronic subprocesses. 
Computer programs modeling these subprocesses, such as the Monte Carlo event generator 
Herwig++ [B + 08], rely on assumptions about the intrinsic motion of quarks inside hadrons, 
see, e.g., Ref. [GSS08]. On the lattice, we can test these assumptions, e.g., how well the 
factorization "/i(ic, k\) oc f\(x)f[ (fejj" is fulfilled or whether a Gaussian function is an 
adequate description of the fcj_-dependence. 

Note, however, that most of the results we present have been obtained with a simplified 
operator, which is particularly easy to realize on the lattice. There is still much debate 
about the precise operator needed to define TMD PDFs that are suitable in the description 
of scattering processes. For a quantitative comparison of results from the lattice and from 
scattering experiments, we will have to go beyond the simplified operator. We shall briefly 
show first steps in that direction towards the end of this work. 



TMD PDFs and GPDs are complementary in the sense that there is no simple transformation among 
them [Ji04]. However, there exist certain non-trivial relations between GPDs and TMD PDFs [Bur04, BH04, 
MMG07]. In principle, it is conceivable to study and fcx-dependence simultaneously in the context of 
Wigner distributions [Wig32, Ji03, BJY04, Ji04]. 



1.1 Outline 
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Numerical lattice calculations are carried out in different stages. The production stage requires 
months on a supercomputer. These expensive computations are justified by the wide range of 
questions that can be addressed with the resulting data. In our case, we profit from the fact 
that the essential building blocks needed to calculate GPDs can be reused to analyze TMD 
PDFs, at a comparatively small additional computational cost. 

This research project is based on configurations generated by the MILC collaboration [B + 01] 
and propagators generated by the Lattice Hadron Physics Collaboration (LHPC) [H + 08]. The 
production of these large files has required an immense computational effort. To store the files 
locally, we have equipped one of our computers with 3.5 terabytes of space on hard drives. 
For the primary steps of our analysis of TMD PDFs, we have developed programs in C++ 
using the Chroma and QDP++ libraries [EJ05] for parallelized lattice computations. To run 
them, we have combined personal computers of our theory department to small clusters with 
the MPICH2 implementation of the Message Passing Interface. For the final analysis we 
have used Mathematica. 



1.1 Outline 



The next chapter will familiarize the reader with the concept of TMD PDFs in the context of 
scattering experiments and factorization theorems. It will also become clear that the correct 
definition and regularization of TMD PDFs is still under debate. 

The third chapter will start with a brief introduction to lattice QCD, before we will specialize 
to the calculation of nucleon structure. We will give some specifications of the ensembles and 
propagators that were kindly provided to us by the MILC and LHPC collaborations, and will 
explain our techniques to estimate statistical uncertainties. 

In the main part of this thesis, chapter 4, we will work with a simplified definition of TMD 
PDFs, suitable for first explorations on the lattice. The TMD PDFs we thus obtain are not 
strictly identical to those used in the literature and for the description of, e.g., semi-inclusive 
deeply inelastic scattering. Nevertheless, they do provide qualitative insights into, e.g., the 
spin structure of the nucleon. A large part of this chapter will be devoted to renormalization, 
a necessary step in order to be able to present results which are independent of the details of 
our lattice calculation. 

The fifth chapter explores future concepts. In particular, we will discuss the results of a test 
calculation which goes beyond our simplified definition of TMD PDFs. 

We summarize in chapter 6. Notes on conventions and some details are provided in the 
appendix. It is planned to create a web page with documentation of the file formats and 
software tools developed for this research project [Mus]. 
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1.2 Basics of QCD 

In general, electroweak interactions have little influence on the structure of hadrons. Thus 
the theory we need for the description of the structure of the nucleon is QCD, which describes 
just quark and gluon fields. For details on QCD, we refer to textbooks such as Ref. [Mut87]. 
The QCD Lagrangian reads 

C QC v[q,q,A)(x) = J2 q(x)(0-m q )q(x)-^F^ a (x)F^(x) , (1.1) 

q=u,d,s,... 

where v are Lorentz indices and where the index a = l..(N^ — 1) refers to the adjoint repre- 
sentation of the color group SU(iV c ), N c = 3. The quark fields u(x), d(x), . . . and u(x), d(x), . . . 
implicitly carry a color index i = l..N c and a Dirac index a = 1..4. In functionals and func- 
tional integrals, we will use the symbols q and q to refer collectively to all quark degrees 
of freedom. The quark masses m u , m^, m s , . . . are fundamental constants. The covariant 
derivative operator 

D lx = d lx -igA li (1.2) 

introduces an interaction of quarks and gluons. Here the gluon field A^ix) is a 3 x 3 color 
matrix, which can be expressed in terms of real fields A^ a (x) according to A^x) = A^ a (x) T a , 
where the T a are the 8 generators of SU(3). 2 The coupling strength is given by the constant 
g. The field strength tensor F^^x) is defined in terms of the gluon field as 

F fiv a = dy,Ai/ a — d^A^a + g f a t, c A^^ A v c , (1-3) 

where f a b c are the structure constants of SU(3). The Lagrangian £qcd is invariant under 
local gauge transformations of the form 

q(x)^q'(x) = W(x)q(x) , (1.4) 
q(x) -> q'{x) = q(x) W\x) , (1.5) 

A^x) A'^x) = W(x) (a^x) - % - W'\x) (d^Wix)) ) W~\x) , (1.6) 

where the unitary 3x3 matrix VF(x) is an element of SU(3). As we will see in section 3.1.5, 
local gauge transformations of the gluonic degrees of freedom will look much simpler in the 
lattice formulation. 



2 In terms of the Gell-Mann matrices A a , the generators are defined as T a = \ a /2. 



Chapter 2 



Nucleon Structure from Deeply 
Inelastic Scattering Experiments 

This chapter briefly reviews the role TMD PDFs play in our understanding of scattering 
processes. The correlators introduced here will be the starting point for our lattice calculations 
described in the rest of this work. However, as we will explain, specifying the correlators 
precisely is a difficult task and still a matter of ongoing research. 



2.1 Experimental Setups 

To determine the structure of the nucleon experimentally, it is advantageous to probe the 
nucleon with interactions that can be described accurately in perturbation theory. Reactions 
of the nucleon with leptons are therefore of primary interest. The interaction between the 
lepton and the parton is mediated by an electroweak particle, i.e. a photon, a W or a Z 
boson. The following reactions are sensitive to TMD PDFs: 

• A + B^l + l + X, the Drell-Yan process, and 

• I + H — > I 1 + h + X, called 1-particle inclusive or 

semi-inclusive deeply inelastic scattering (SIDIS). 

The two processes given above are illustrated and explained in Fig. 2.1. Note that we re- 
strict ourselves to the leading contributions, where only a single electroweak gauge boson is 
exchanged. 

In the following, let us have a closer look at SIDIS. For reasons of clarity, we restrict ourselves 
to the special case that the leptons are electrons e~, the incoming hadron is a nucleon, and 
the exchanged virtual gauge boson is a photon 7* . 

The kinematics of the reaction is depicted in Fig. 2.2a, see Ref. [BDDM04]. The momenta 
Pi and Py of the in- and outgoing electrons together with the nucleon momentum P define 
the lepton plane. We choose to align the nucleon momentum P and the photon momentum 
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(a) (b) 




Figure 2.1: Two examples of processes sensitive to TMD PDFs. We draw the leading contri- 
butions, in which a single electroweak gauge boson (wiggled lines) is exchanged. 

(a) The Drell-Yan process. Two hadrons A and B collide to form a lepton pair I, I and a 
number of other particles subsumed in X. 

(b) The SIDIS process. A lepton I scatters off a hadron H (typically a proton), leading to 
the production of a hadron h as part of a jet. Apart from h, all particles in this jet and in 
the debris of H are collected in X. 



q = Pi — P[i with the z-axis 1 . Then the momentum of the produced hadron Ph can have 
components Ph± in the xy-plane, transverse to the nucleon momentum. Thus we can define 
an azimuthal angle (ph between Ph± and the lepton plane. For a polarized nucleon target, the 
transverse spin components of the nucleon form an angle (j>s with the lepton plane. 

In the SIDIS cross section 

the lepton tensor L^ u is calculable in perturbation theory. All the non-perturbative informa- 
tion related to hadron structure is encoded in the hadron tensor 

W^(P, q, P h ) = 5^(q + P-P h -P X )Y, (N(P, S)\ J»(p) \X h(P h , S h )) 

x 

x (Xh(P h ,S h )\ J"(0) \N(P,S)} 

-J4 



/ E (N(p,s)\j»(e)\xh(p h ,s h )} 



x 

x (Xh(P h ,S h )\r(0)\N(P,S)) , (2.2) 

where \N(P, S)) is a nucleon state with momentum P and spin S, h(Ph, Sh) is a hadron with 
momentum Ph and spin Sh, and where J^(£) is the electromagnetic current of the quarks at 



A common choice of this type is the Breit frame, with the only non-vanishing component of the photon 
four-momentum q 3 = —Q. 



2.2 A Parton Model inspired Factorization Ansatz 
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(a) (b) 




Figure 2.2: (a) Kinematics of SIDIS following the Trento conventions [BDDM04]. 
(b) Simplified factorized tree level diagram of the hadron tensor in SIDIS. 



position I. Any one of the two matrix elements on the right hand side of eq. (2.2) corresponds 
to the gray blob in Fig. 2.1b. Note that the process is called semi-inclusive because we sum 
over X, which is only part of the final state \X h(Ph, Sh))- 

2.2 A Parton Model inspired Factorization Ansatz 

How can we extract information about the structure of the nucleon from the hadron tensor? 
We need to decompose ( "factorize" ) the reaction further into perturbative (hard) and non- 
perturbative (soft) subprocesses. A simplified sketch of a factorized hadron tensor is given in 
Fig. 2.2b. The left and right halves of the diagram are mirror images. Each half is related 
to one of the matrix elements on the right hand side of eq. (2.2). Let us concentrate on the 
left half of the diagram. The central assumption we want to make is that the virtual photon 
7* couples with the electromagnetic current of a single quark of the nucleon, which initially 
carries momentum k. We want to look at cases where the virtuality Q 2 = —q 2 of the photon 
is large, Q 2 3> m 2 N . In this case the quark-photon-interaction is a hard process, i.e., it can be 
treated perturbatively and it is drawn outside the gray blobs. The struck quark cannot appear 
as a real particle in the detector because it carries a color charge. According to the confining 
property of QCD, color charges cannot be isolated from each other over distances more than 
about a femtometer. Thus the struck quark undergoes a process called hadronization, or 
fragmentation, depicted in the upper gray blob of the diagram and parametrized in terms 
of fragmentation functions (also called decay functions). Here quark-antiquark pairs emerge 
from the vacuum and guarantee that the particles appearing in the detector are color neutral 
objects like hadrons. A jet of particles forms, of which only one hadron is detected individually. 
The rest of the jet (upper blob) and the remnants of the incoming nucleon (lower blob) are 
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subsumed in X and are typically ignored in the analysis of the experiment. 

The proposed decomposition of the hadron tensor in Fig. 2.2b corresponds to an early factor- 
ization ansatz of the cross section suggested by Collins [Col93] 2 : 

d^P h = JfW £ Jd 2 k x fi, q (x B ,kr,Q 2 )^D h>q (z,k x + q± ;Q 2 ). (2-3) 

1 n q=u,d,s,... 

Here we have used the following kinematic variables: 

In eq. (2.3), the cross section a collects the hard pieces of the reaction, i.e., the short-distance 
part of electron-quark scattering. Dh jQ (z, k± + q ± , Q 2 ) is a fragmentation function for quark 
type q and corresponds to the upper blob in Fig. 2.2b. Finally, fi <q (xB, k±; Q 2 ) describes the 
density of quarks q in the nucleon. It corresponds to the lower blob in Fig. 2.2b and is the 
prototype of a TMD PDF. 

Why do we parametrize the quark distribution fi <q in terms of xb and k±? Let us understand 
this in the parton model, where we assume that the moment of interaction with the lepton 
is so short that the nucleon can be described as a collection of free particles, the "partons" 
(quarks or gluons). Consider, for example, the rest frame of the incident lepton, Pi = 0. In 
this frame, the momentum of the nucleon will be very large, which we can write as P + 3> mjv 
using the light cone coordinates specified in section A. 2. The entire momentum transfer q is 
passed on to a single parton, a quark with initial momentum k. Because we work in a large 
momentum frame of the nucleon, the relative magnitudes of the components of k will behave 
as k + : k± : k~ ~ P + /mjy : 1 : rriN/P + . We neglect the suppressed momentum component 
k~ in the kinematics of the process. To specify k + , we introduce the dimensionless variable 
x = k + /P + . Now, let us require that the parton be not too far off-shell, i.e., k 2 ~ and 
(k + q) 2 rj 0. This means ~ k 2 + 2k-q + q 2 ~ 2xP-q — Q 2 , from which follows x ~ xb- Thus 
the so-called Bjorken scaling variable xb can be approximately identified with the longitudinal 
momentum fraction x of the parton in the nucleon. 

We have already remarked that k± is suppressed by a factor mjv/P + as compared to the 
longitudinal quark momentum. In many cases, the average intrinsic transverse momenta of 
quarks in the nucleon and the final state are negligible in the kinematic description of a 
scattering process. In such cases, one can work with the conventional, "integrated" PDFs 
and fragmentation functions like fi, q (x;Q 2 ) and Dh jq (z; Q 2 ), compare section 2.3.2 below. 
However, when we consider SIDIS for a small transverse jet momentum (|<jf_J *C Q), we 
can no longer ignore intrinsic transverse momenta, and we need "unintegrated" , transverse 
momentum dependent distributions. Just like an integrated parton distribution function 
fl >q (xB] Q 2 ), the transverse momentum dependent distribution fi <q (x B , k±; Q 2 ) follows an 
evolution equation in Q 2 , which has been studied in Ref. [CT06]. 



2 To keep the discussion simple, we quote the unpolarized cross section. 
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To promote the factorization ansatz to a factorization theorem, one can stick to the following 
strategy [CSS88]: First, one evaluates the cross section in perturbation theory. To make this 
possible, incoming and outgoing hadrons are replaced by quarks. Then, one shows that the 
factorization ansatz is indeed valid for the perturbative calculation. This involves proving 
that the leading contributions come from Feynman graphs with momenta that are either 
far off-shell, or inside non-overlapping "leading regions" , each of which is attributed to one 
of the non-perturbative contributions f q and D^ q. One can now identify the contributions 
from f q and P>h,q in the cross section and divide them out. For the remaining short-distance 
contribution a it is irrelevant whether we are using quarks or hadrons as external particles. 
Thus finally factorization is established and a perturbative expression for the short-distance 
part <7 is available. 

However, due to the important role of gluons, it turns out that Fig. 2.2b and eq. (2.3) need 
to undergo modifications before they can be used as an ansatz for a factorization theorem. 
We shall discuss these difficult issues later. 



2.3 The Fundamental Quark-Quark Correlator 
2.3.1 Definition 

Let us stick to the naive picture of the previous section for a moment and look once more 
at the lower blob in Fig. 2.2b. It is a quark-quark correlation function of the nucleon. In its 
general form, we write it as 



$f P \k, P,S)= ™ e -*4 (iV (P, S)\ q(£) T°p Ufa) q(0) \N(P, S)) , (2.5) 



where the r op is a Dirac matrix. The Wilson line U[Cf\ as defined in eq. (A.l) runs along a 
path Ci from I to the origin. It is necessary to connect the quark operators via a Wilson line 
to ensure gauge invariance. The meaning of the Wilson line, also called gauge link, and the 
choice of the contour Ci will be discussed in the sections to follow. 

We have already mentioned that k~ is negligible in the large momentum frame of the nucleon. 
Therefore, we are going to average over it in the correlator. We can also understand this as 
follows: Imagine the reaction takes place in the xy-plane. The time it takes for the nucleon 
to traverse this plane is proportional to m,]\r/P + . Boosting to the nucleon rest frame, we find 
that the region of spacetime locations I where the reaction takes place has a negligible extent 
in the "+"-direction, suppressed by [m^ / 'P + ) 2 . Therefore, the important information for the 

~ [Pop] 

description of the high momentum collision is given by the correlator % 1 {I, P, S) at l + = 0. 
After the Fourier transform, this is equivalent to an integration over k~ . Thus we define 





(2.6) 
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with the longitudinal momentum fraction x = k + /P + as introduced before. The calculation 

rpopi 

of the transverse momentum dependent correlators & q (x, k±; P, S) will be the focus of the 
study at hand. In order to state the results in a form independent of an explicit choice of a 
frame of reference (in particular, independent of P + ), these correlators are parametrized in 
terms of TMD PDFs, see section 2.6. TMD PDFs are profile functions that kinematically 
only depend on x and k\. One TMD PDF we have already met is fi q (x, k±). It appears in 

(x, k ± ; P, S) = fi }q (x, k\) + (spin dependent terms) (2.7) 

and is commonly interpreted as the density of quarks of flavor q in the nucleon, see section 2.5. 
It is easy to understand why the correlator with r op = 7+ is so important: As in the case of the 
quark momentum k, the "+"-component of the current (N(P, S)\ q^q \N(P, S)} is enhanced 
by a factor P + /m^- Note that & q and fx q are also functions of renormalization and 
factorization scales fi, Q 2 , . . ., as we will see later. For the moment, we omit these additional 
variables in our notation. 



2.3.2 Relation to PDFs 

If we integrate not only over k~ , but also over k± in eq. (2.5), the components l + and £± 
become zero: 



d>P(x;P,5) = j dk~ jd 2 k ± &P(k,P,S) 



^e~ ix P+£ ~ l(N(P,S)\ r«-)r op W[C„_] g (0) WP,S)) . (2.8) 

Z7T Z 



This correlator is used to introduce conventional, fcj_-integrated PDFs such as fi,q(x) 

4 7+1 (^ 

such as 



$f 4j (x;P,5). Thus, naively, one would think that relations between PDFs and TMD PDFs 



" fi, q (x) = J d 2 k ± /i, 9 (z, fci) " (2.9) 

hold. However, the integral above is undefined. This is due to the behavior of fi,q(x, fcjj at 
large k±, where perturbation theory is applicable. As explained in, e.g., in Ref. [BBDM08, 
Die08], one finds: 



1 



fl,q{x, k±) ~ T^h,qi x ) for k ± -> 00 • ( 2 - 10 ) 

k ± 

The interesting, non-perturbative information is encoded in TMD PDFs at low k±. Intro- 
ducing a cutoff in eq. (2.9), 

Q) = I d 2 k ± f 1<q (x, k\) (2.11) 

J\k ± \<Q 

regularizes the integral and leads to a Q-dependence that corresponds to the DGLAP evolution 
[GL72, Par 74, Dok77, AP77] of conventional, integrated PDFs. However, for it is 

no longer evident that the correlator simplifies as in eq. (2.8), see also the discussion in [Col03]. 



2.3 The Fundamental Quark-Quark Correlator 



17 



The approaches mentioned in sections 2.4.4.2, 2.4.4.3, and 2.4.4.4 promise to elucidate the 
relation between PDFs and TMD PDFs. For the calculation of integrated PDFs in practice, 
the divergence in k± does not pose a problem, because the correlator eq. (2.8) is renormalized 
and evaluated directly, without the detour via TMD PDFs. 



2.3.3 Mellin Moments 

Just as for PDFs, it can be useful to analyze the x-dependence of the correlator in terms of 
so-called Mellin moments. The n th Mellin moment is defined as 

$pW(fe i; ?,5) = J' dx x n - 1 *P(x,fe ± ;P,5) 

= Jdxx n ~ 1 ($f°^{x,k ± ;P,S) + {-l) n - 1 $f° P] {-x,k ± -,P,S)} . (2.12) 

[pop] 

Given that & q (x,k±; P,S) vanishes for \x\ > 1, the integration limits in the definition 
above are inessential. As will be explained in section 2.5.1, we associate negative values of x 
to antiquarks. Thus the second line of the above equation indicates that the Mellin moments 
introduced this way are a combination of quark and antiquark distributions. Of particular 
interest to us is the first Mellin moment, where we simply integrate over x = k + /P + : 

= ^J^^^r ] (i,P,S)\ e+ ^ . (2.13) 

The quark separations £ appearing in the correlator now lie in the transverse plane and 
are purely spatial - a good premise for lattice calculations. The correlator describes the 
distribution of quarks (and antiquarks) in transverse momentum space, irrespective of their 
longitudinal momentum. For example, using the relation fi,q(x, k\) = -f ljQ (-x, k\) [TM95], 
we have in the unpolarized case 

^ +](1) (k±;P,S) = fi}{k\) + (spin dependent terms) , 

fxl&l) = I dx f hq (x, k 2 ± ) - f 1 dx f 1>q (x, k 2 ± ) . (2.14) 
Jo Jo 

Thus we interpret as the difference of two fcj_-dependent densities, namely the quark 

density dx fi, q {x,kl) and the antiquark density Jq dx fi jq -(x,kl). We present lattice 
results for the first Mellin moment in section 4.6. 
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(a) (b) 

gov cov+£ 




Figure 2.3: (a) Cut SIDIS diagram with gluons at leading order in M/Q. Longitudinally 
polarized gluons are indicated by "+", transversely polarized gluons by "_L". There can 
be any number of longitudinally and transversally polarized gluons. Note that the total 
momentum carried by the quark and the gluons leaving the blob is k. 

(b) Diagram of the quark-quark correlator with gauge links and gluons. In SIDIS, the link 
direction v is either identical or close to h-, depending on the regularization prescription. 
Partial cancellation of the sections of the gauge link which run to transverse infinity has been 
disputed in Refs. [CS08a, CS08b]. 

2.4 Improved Definitions of TMD PDFs 

2.4.1 A Starting Point: The Straight Wilson Line 

An naive guess for the Wilson line U[Ce] in eq. (2.5) is a straight line from £ to 0, i.e. U[Ci\ = 
U[£, 0] in the notation of eq. (A. 2). The resulting correlator is gauge invariant, and serves us 
for first exploratory studies on the lattice, see chapter 4. However, definitions of TMD PDFs 
designed for the description of real scattering experiments require more complicated gauge 
link structures, see below. 

2.4.2 The Physical Role of the Wilson Line 

The Wilson line U[Ce] in the quark-quark correlator arises naturally from diagrams involving 
gluon loops that contribute at leading order in mj^/Q. In the following, we briefly moti- 
vate this statement in the context of SIDIS. Details can be found, e.g., in Ref. [Pij06] . We 
will find that the Wilson line represents "final state interactions", namely soft gluons ex- 
changed between the remnants of the nucleon and the hadronizing parton (compare, e.g., 
Refs. [BHS02a, BHM+02]). 

In the attempt of proving a factorization theorem for SIDIS, diagrams of the type shown in 
Fig. 2.3a turn out to be important. Consider j gluons with momenta pi, p%, . . ., pj. From 



2.4 Improved Definitions of TMD PDFs 



19 



the quark propagator before the i th gluon vertex, we get a denominator of the form 

1 1 1 



(k + q-pi -pi+i - ... ~Pj) 2 -m 2 + ie -2q n_ • (p; + . . . + pj) 



le 



(2.15) 



The approximation above is called eikonal approximation (see Ref. [CS81] for details) and is 
applicable for soft gluons (the pi are small) and quarks almost on-shell ((k + q) 2 — m 2 ~ 0). 
Now let us look at a straight Wilson line running from the origin to infinity along a four- vector 
v: 

U[oov, 0}=V exp (-i g J dX v-A(Xv)j = V exp g J dX v-A(Xv)"j 



x f'OO 

i + Y,^) j / dA i 



dXj v-A(vXj) ■ ■ ■ v-A{vX 



i + V I d Pl f d Pj gv ' A< ^ p ^ gy-Api) (2 16) 

fiJ (27T) 4 '"; (2vr) 4 vpj-ie " ' v ■ ( Pl + . . . + Pj ) - ie ' ' 



Here V denotes reverse path ordering and A(p) = f d 4 x e tpx A(x). The last line has an 
interpretation in terms of Feynman rules [CSS88, CS82, CFP80, BFK80]: The denominators 
{v-(pi + . . . + Pj) — ie) -1 are displayed as double lines and are called eikonal lines. The gluon 
vertex is proportional to v^, i.e., only gluons polarized along v can couple to the eikonal lines. 
The denominators in the last line of the equation above remind us of eq. (2.15) if we set 
v = n_. Indeed, we are able to encode the exchange of the longitudinally polarized gluons 
(coupling with A + and marked with a "+" in Fig. 2.3a) in our TMD PDFs by including a 
Wilson line Li[oov 1 0] in the definition of our quark-quark correlator P, S). 

There is a trick frequently used in the literature to hide the Wilson line. We can define quark 
fields with a "gauge history" 

^ q:V {i)=U[oov + £,i\q{t) (2.17) 
and write our quark-quark correlator in eq. (2.5) as 

&^(£,P,S) = \ (N(P,S)\ VvWr op <M0) \N(P,S)) . (2.18) 

However, for v = n_ this trick only works in gauges where the transverse gauge fields 
^4_i_(oon_ + i±) vanish at light cone infinity [JMY05]. In other gauges (including light- 
cone gauge A + = 0), it has been noticed that there are also leading contributions from 
transversely polarized gluons with negligible momentum in "+"-direction [BJY03, BMP03, 
Pij06]. We can absorb these effects into the soft correlator using the transverse Wilson line 
U[ooh- + oo£±, oon_]. Figure 2.3 shows the resulting quark-quark nucleon correlator for 
SIDIS. Together with the other side of the cut diagram, the complete Wilson line in eq. (2.5) 
reads 

U[C e ]=U[£, oov + £, oow + oo^] ^[oov + ooIl, oov, 0] , (2.19) 

where v = n_ and where we have used the notation introduced in eqns. (A. 2) and (A. 3). 
In some works, it has been assumed that the transverse sections of the Wilson line cancel 
partially, so that 

U[C t ] =U[£, oov + £ ± , oov, 0] . (2.20) 
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The shape of this Wilson line is a staple extending out to infinity. Note, however, that 
Refs. [CS08a, CS08b] claim that the cancellation of transverse sections is incorrect due to 
additional divergences produced by a cusp in the Wilson line at transverse infinity (see sec- 
tion 4.4 for the discussion of renormalization properties of Wilson lines). In any case, the 
Wilson line now forms a continuous connection between the two quark fields, so that we end 
up with a gauge invariant definition of TMD PDFs. 

In our considerations above we have focused on the SIDIS process, where we obtain a Wilson 
line that corresponds to final state interactions. It should be remarked that for the Drell-Yan 
process, the lightlike part of the Wilson line runs in the other direction as compared to SIDIS, 
i.e., v = —h- [Col02]. In this case, the Wilson line represents "initial state interactions" 
(compare, e.g., Ref. [BHS02b]). In general, an incoming quark turns into a Wilson line 
coming in from — oon_, while an outgoing quark creates a Wilson line out to oon_ [BMP06, 
Pij06, CKKL05]. 

2.4.3 Rapidity Divergences in the Lightlike Wilson Line 

In the eikonal approximation eq. (2.15), the ejected quark is treated like a massless particle 
moving along the "—"-direction. It was realized [CS82] that this procedure removes a physical, 
process dependent cutoff, and thus leads to a severe divergence, sometimes termed rapidity 
divergence, see, e.g., Ref. [CRS08]. The cause of the divergence can be traced back to gluons 
with unphysically large momentum in the "—"-direction [CS82]. The divergence cannot be 
regularized by the introduction of a gluon mass or with dimensional regularization. Two types 
of strategies have been developed to handle this divergence: 

1. Following a suggestion by Collins, Soper and Sterman [CS81, CS82] in the context of 
fragmentation functions, the longitudinal gauge link can be placed slightly off the light 
cone, i.e., with v not exactly equal to n_. The distributions defined in this way depend 
on Q = (P • v) 2 /\v 2 \ ~ (P + ) 2 \v~ /2v + \, which acts as a cutoff parameter and is chosen 
large but finite. An evolution equation can be derived to describe the dependence of 
the distributions on £, see e.g. eqns. (6.4) and (6.6) in Ref. [CS81]. 

2. The divergence can be cancelled by a subtraction factor, introduced as vacuum expec- 
tation values of Wilson lines. In general, in such an approach the parameter controlling 
the rapidity cutoff will appear in the subtraction factor. 

A review of these approaches can be found in Ref. [Col03] . 

2.4.4 Definitions of TMD PDFs Proposed in the Literature 

Let us look at some concrete proposals in the literature how the Wilson lines can be arranged 
to obtain well-behaved definitions of TMD PDFs. 
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2.4.4.1 The Factorization Formula by Ji, Ma and Yuan 

Ji, Ma and Yuan, [JMY05, JMY04] have proposed a factorization prescription where both 
Wilson lines slightly off the light cone and subtraction factors appear. Let us study their 
results for SIDIS. In their definition of TMD PDFs, which we illustrate in Fig. 2.4, they 
introduce two slightly timelike link directions v = v + h++v~h- ~ n_ and v = v + h++v~h- ~ 
n + . To factorize the hadron tensor, Ji, Ma and Yuan take the most general setup of leading 
regions into account, as depicted in Fig. 2.5. As an example, consider the unpolarized leading 
structure function F in the hadron tensor W^ v = — ^g^F(xB, z, Ph±, Q 2 ) + • • •• In factorized 
form, the Fourier transform of F with respect to Ph± reads 

F(x B ,z,£ ± ,Q 2 ) = ^2 e gh,q( x B,z£±;^ 2 ,Q 2 p,p)Dh : q(z,e±;fi 2 ,Q 2 p,p) 

q=u,d,s,... 

x S(£±;p 2 ,p) H(Q 2 ,p 2 ,p) . (2.21) 

Here e q is the electric charge of quark q, and p is a renormalization and factorization scale. 
The dependence on which parametrizes the direction of v, is hidden in the dependence on 
Q 2 and the parameter p = v~v + /v + v~ via a special choice of coordinates. The soft factor 
S appears not only in the equation above, but also as subtraction factor in the definition of 
the fragmentation function Dh q and in the definition of the TMD PDF: 



f hq (x,£ ± ;p 2 ,Q 2 p,p) = 



S(£±,p 2 ,p)J (2tt) 



tP+l- 



x^(N{P,S)\ q(£)U[l,oov + l]j + U[oov,0]q(0) \N(P,S)) ^ . (2.22) 

Since Ji, Ma and Yuan work in a gauge where gauge fields vanish at infinity, the transverse 
sections of Wilson lines at infinity need not be specified. Thus up to the soft factor, this 
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Figure 2.5: The leading regions for SIDIS after soft and collinear factorizations. 




Figure 2.6: Arrangement of Wilson lines in the subtraction factor appearing in the definition 
of TMD PDFs by Hautmann [Hau07]. 



definition is equivalent to the one in eqns. (2.5), (2.7). The soft factor reads 



S(£±;v 2 ,p) = (0|tr U[-oov + £±, £±, oov + £±] U[oov, 0, -oov] |0) . (2.23) 



Ji, Ma and Yuan provide evolution equations for the p- and ^-dependence (see also Ref. 
[IJMY04]), and give arguments that their formula is valid to all orders in perturbation theory. 
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2.4.4.2 Subtracted TMD PDFs by Hautmann, Collins and Metz 

Hautmann [Hau07] follows suggestions put forward in Refs. [CHOO, CHOI, CM04] and proposes 
a quark distribution of the form 



fi, q (x,£±;Cu) 



dt 



cP+£~ 



(27T) 

N(P 

(0| tr U[oou + ra„£~, h-l~, oon_] U[ooh-, 0, oou] |0) 



x-(N(P,S)\ q(£)U[£,ooh- + £} 7 + U[oon^,0}q{0) \N(P,S)) 



(0| tr U[oou + £, £, oon_ + £] U[ooh_, 0, oou] |0) 



£+=0 



(2.24) 



Again, a gauge with vanishing gauge fields at infinity is employed. The subtraction factor in 
the last line of the above equation is illustrated in fig. 2.6. Its purpose is to cancel "endpoint 
singularities" occurring in TMD PDFs for x — * 1. In contrast to the approach by Ji, Ma 
and Yuan, the vector v employed here is exactly n_. The auxiliary non-lightlike direction 
u = u + h + + u~n_ appearing in the subtraction factor introduces a dependence on the reg- 
ularization parameter Q u = (P + ) 2 u~ /2u + . However, the subtraction factor cancels and the 
dependence on £ u disappears for the fcj_-integrated parton distribution function, which one 
obtains by setting £± = in eq. (2.24). Thus the definition by Hautmann promises to provide 
a relation between TMD PDFs and regular, fcj_-integrated PDFs. 

2.4.4.3 TMD PDFs in Soft-Collinear Effective Theory by Chay 

Another proposal [Cha07] promising a relation to regular PDFs employs soft-collinear effective 
theory (SCET) to regularize the parton distributions, using exclusively light-like Wilson lines. 
The definition of the quark distribution contains the subtraction factor depicted in fig. 2.7 
and reads 

h, q {x,k ± , K ) -J du J (2vr)2(27r)2 e 

x (N(P, 5)| Xn + ^±) 5 ^ - ^ + )^Xn + (0) \N(P, S)) 

1 (0\trU[ooh + +£, £, -oon_ + £] U[-ooh-, 0, -ooh+) |0) . (2.25) 



N c 



e+=o 



Here Xh+ are quark fields within SCET. Their definition includes a "gauge history" analogous 
to eq. (2.17). The appearance of u and the large momentum label operator is a SCET 
specific feature. The additional scale k in the above equation is related to soft gluon emission. 
Note that the Wilson line in the soft factor starts at — oofi+ but ends at con+ +£. Thus there 
is no obvious way to close the Wilson line to a manifestly gauge invariant loop. 



2.4.4.4 Definition of TMD PDFs by Cherednikov and Stefanis 

As already mentioned, the authors of Refs. [CS08a, CS08b] deduce from an analysis of anoma- 
lous dimensions that the Wilson lines in eqns. (2.19) and (2.20) are not equivalent. They 
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Figure 2.7: Arrangement of Wilson lines in the subtraction factor of Chay's definition of TMD 
PDFs in SCET [Cha07]. 

suggest a definition of parton densities with lightlike Wilson lines (v = n_) and include soft 
subtraction factors with cusps at and £: 



Here the mass scale r\ is hidden in the regularization of rapidity divergences with a pole 
prescription and plays a role similar as £ in section 2.4.3. The renormalization scale [i is 
needed for dimensional regularization. The choice of the transverse direction u± is completely 
arbitrary. An illustration of the Wilson lines is given in fig. 2.8. As in Chay's approach, there 
is no obvious way to complete the Wilson lines such that the soft factors become manifestly 
gauge invariant expressions. Cherednikov and Stefanis show that their TMD PDFs fulfill - 
at least formally - the simple relation to integrated distributions eq. (2.9). 

2.4.5 Wilson Line Self-Energy 

As will be discussed in section 4.4, the Wilson line exhibits a divergence dependent on its 
length L due to self-energy graphs. It can be removed by a factor exp(— 5m L), where 8m is 
a renormalization constant which vanishes for dimensional regularization, but not for cutoff 
schemes as in lattice QCD. The divergence may be eliminated by a subtraction factor included 
in the definition of the TMD PDF. In accordance with the suggestion in Ref. [Col08] , consider 




x - (N(P, S)\ q{£) U[£, oon_ + £ ± , oon_ + oou ± + £ ± ] 7 4 

xU[oon^ + oou±, oon_, 0] q(0) \N(P,S)) 

x — (0| tr U[— oon+, 0, oon_, oon_ + oo£±] |0) 

x — (0| tr W[oon_ + oou± + £, oon_ + £, £, ooh + + £] |0) 



(2.26) 
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Figure 2.8: Arrangement of Wilson lines in the definition of TMD PDFs by Cherednikov and 
Stefanis [CS08a, CS08b]. The cross at the end of the transverse link indicates a non-smooth 
connection of the line at infinity. 



the definition 

f 1 , q (x,kr,yn)= ^L e - P+e -+^ Inn 



(2vr) J (2vr) 2 v^°° 
1(N(P,S)\ q(£)l + W, r)v + £, rjv, 0] q(0) \N(P,S)) 



X 2 yjtoop^l) 



(2.27) 

£+=0 



The numerator in the second line is our $t 7+ ], however the Wilson line U\Ci\ does not run all 
the way out to infinity. The rapidity parameter y n symbolically reminds us that we must take 
care of rapidity divergences, e.g., by choosing a non-lightlike v. The Wilson loop (oop(L,£) is 
given by: 

tbop(ri,£) = — (0|tr U[rjv, -r)v, -rjv + £, rjv + £, rjv)\0) . (2.28) 

This loop is twice as long as the link path in the numerator of the second line of eq. (2.27). 
Thus the square root in the denominator cancels the self-energy divergence. Note that the 
loop expectation value becomes //-independent for v = n_. 

We will try out the idea presented above on the lattice in section 5.2.2, with a setup of link 
paths as illustrated in Fig. 5.3. 

It might be interesting to generalize the idea of the loop subtraction factor to cusped loops 
of the form 

Coop(rj, £) = ^v- (0| tr U[0, -qv, rjv + £, £, rjv + £, ryu, 0] |0) , (2.29) 

where we have introduced an auxiliary direction v. The loop has now approximately the 
shape of a folded rectangle, see Fig. 2.9. With an appropriate choice of v, this loop structure 
cancels the self-energy divergence and, taking into account that transverse pieces can be left 
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Figure 2.9: Cancellation of self-energy contributions in the definition of TMD PDFs, as 
suggested in Ref. [C0IO8] but generalized to two different directions v and v in the Wilson 
loop subtraction factor. 



out in the limit 77 — > 00 in appropriate gauges, shows similarity to the subtraction factors 
introduced by Ji, Ma, and Yuan (cf. section 2.4.4.1) as well as Hautmann, Collins and Metz 
(cf. section 2.4.4.2). However, these authors do not put their subtraction factors under a 
square root. 

As a side remark, note that the middle section U[rjv + £, r]v] of the gauge link in eq. (2.27) is 
not exactly transverse if £~ 7^ 0. We could have maintained a transverse middle section by 
taking the gauge link U[£, rjv + £±, r]v, 0] instead. In the limit rj — » 00, the two versions are 
formally equivalent, because oov + £ = oov + h + + (oov~ + £~)h- + £± = cov + £±. If we were 
to insist on the gauge link U[£, r]v+£±, r]v, 0], an appropriate Wilson loop subtraction factor 
(oop(r],£) would look slightly more complicated, and the simple transformation properties of 
the gauge link under discrete symmetries in eq. (2.48) below would not hold exactly any more. 



2.4.6 Remarks 



Quoting from a recent proceedings article by Collins [C0IO8] - "to allow non-perturbative 
methods in QCD to be used to estimate parton densities, operator definitions of parton 
densities are needed that can be taken literally". As we have seen, there is now a variety 
of proposals available for transverse momentum dependent parton correlators. We hope that 
from these ideas a generally accepted definition of TMD PDFs will emerge that remains 
finite-valued and well-defined beyond perturbation theory. 
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2.5 Probability Interpretation of TMD PDFs 

2.5.1 The Parton Picture using Light-Front Quantization 

We usually think of TMD PDFs as probability distributions of quarks inside the nucleon. This 
notion has its origin in light-front quantization [KS70, Dir49, BPP98, Bur96]: As mentioned 
in section 2.3, the reaction at high momentum is sensitive to the nucleon structure close to 
the plane £ + = 0. The Dirac spinor q(£) on that plane has only two independent dynamical 
components. These "good components" are projected out according to q^(£) = ^7 _ 7 + l(£)- 
The helicity projected Fourier transforms 

q {+] , X (k\k ± ) ^ ^(1 + A 7 5 ) e^-'^x g(+) ( rL+4) (2.30) 

are written in terms of quark annihilation operators b\^ q {k + ,k±) and antiquark creation op- 
erators q (k + , k±) of helicity A and quark flavor q according to 

where f7( + )^(A; + , k±) and V( + ^\(k + , k±) are appropriate Dirac spinors, and where we attribute 
anticommutation relations to the creation and annihilation operators as in Ref. [Bur96]. In 
light-front quantization, the distinction between matter and antimatter is made according to 
the sign of k + (or x = k + /P + ). Note that the sign of k + is equal to the sign of k° if the 
dispersion relation k 2 = m q is satisfied. It is now easy to show that for x > 

V b{ q (xP+,k ± )bgxP+,k ± ) \N(P,S)) _ (2vr)3 [7+] (r , pq) . . 

2^ (n(p,s)\n(p,s)) --pr*u=iK x > k ±>r,b) , V-u) 

A 

i.e., we obtain our usual quark correlator ^ 1+ \x, k±; P, S), though we have to set the gauge 
link U[Cf\ to unity. The operator 6^ q (xP + , kj_) b\^ q (xP + , k±) simply gives the number density 
of quarks with helicity A, momentum fraction x and transverse momentum k±. Thus the 
corresponding TMD PDF fi t q(x,k±) has an interpretation as a density of quarks 3 . 

In order to ensure that the gauge link on the right hand side of eq. (2.32) is unity, we have to fix 
an appropriate gauge. The most natural gauge to use together with light-front quantization 
is the light cone gauge, A + = 0, which, however, is affected by the same divergences as 
lightlike Wilson lines, see section 2.4.3 and Ref. [Col03]. We caution the reader that there are 
conceptual difficulties related to the application of light-front quantization to QCD, see, e.g., 
Ref. [Bur97, Bur96]. 



3 The prefactor (2tt) 3 /P+ is a consequence of the convention to operate on our distributions with integrals 
J dx and J d 2 k± instead of the usual momentum integrals J dk + /(2ty) and J d 2 k j_/(2n) 2 . 
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2.5.2 Positivity of the Quark-Quark Correlator 



Ref. [BBHMOO] gives very useful bounds on TMD PDFs. In particular, it is shown that 
fl,q(x,k±) is positive. Here we reexamine the positivity argument, paying special attention 
to the role of the Wilson line. 



2.5.2.1 Concept of the Proof 

Consider Dirac matrices r op for which a matrix T d exists such that 

l 7 °r op = (r d ) t r d . (2.33) 

For example, to analyze fi >q (x, £Cj_), we select r op = 7 + , then T d = 2 _1 / 4 (^7 _ 7 + ). In the 
following, whenever we work with a specific gauge fixing condition G, we will add an index G 
to the state vectors. In all our considerations, we restrict ourselves to the plane £ + = and to 
x > 0. The Wilson line running from I to the origin in the quark-quark correlator eq. (2.5) will 
have to be split into two pieces according to U[C(\ = U[Ci/]U[C 2t i] with suitable definitions 
of the paths C\£ and C 2 / discussed in the following sections. Let {|n)} be a complete set 
of momentum eigenstates with momenta p n , normalized according to (n\n') = 5 nri i. The 
quark-quark correlator can now be written as: 

i>[ r ° p ](^) = {N{_P,S)\q\£)U[C XJl } (r d )t T d U[C 2t e]q(0)\N(P,S)) = 

G (iv(p,5)|gt WW [ Cv ](r d )t (X>>g cwl r d u[c 2/ ] q (o)\N(p,s)) G = 

( G {n\ T A U[C^ q{i) \N(P, S)) G )* ( G (n| T d U[C 2/ ] #) \N(P, S)) G [ 

n 

E ^ (P - p " } - £ ( G (n| T d W[C M - epq(0) \N(P, S)) G ) * ( G (n| r d ^[C 2/ ] 5(0) \N(P, S)) G ) . 

n 

(2.34) 

In the last step, it is required that the gauge fixing condition G is translation invariant. 
Otherwise G would change into a different condition G' in the matrix element on the left. 
Suppose that our prescription to split the Wilson line fulfills 

U[C U - l] ] = U[C 2>e ] = U[C 3 ] , (2.35) 
where the path C3 is ^-independent. Then 

$[r-] (€) = ^ e i(P- P «H £| gH r%U[C 3 ] ljqf3j (0) \N(P,S)) G 2 . (2.36) 
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For clarity, we have made color indices i,j and Dirac indices a, (3 explicit in the last line. If 
the condition G is independent of £, we can now Fourier-transform to momentum space 



<fl r ° p ](x,/t ± ) 



dt 



d 2 £ l 



2vr J (2tt) 2 
dt_ r d 2 £ 
J W) 



e+=o 



L e i(P-Pn-k)-£ 



G (n\T d U[C 3 ]q(0)\N(P,S)) G 



e+=o 



£ IgH r d W[C 3 ]g(0) \N(P,S)) G 6(p+-(l-x)P+)6 ( - 2 \p n± + k x ) >0. (2.37) 



Apart from proving positivity, the last line gives another inspiration for an interpretation as 
a parton distribution. The delta functions restrict the momenta of the states \n). Obviously 
these states carry the longitudinal momentum of the nucleon, except for a momentum fraction 
x, which is missing, and a transverse momentum —k±. This indicates that the quark field 
operator T d U[Cs]q(0) annihilates a "parton" with longitudinal momentum xP + and transverse 
momentum k±. 

Ref. [BBHMOO] uses the argument above for various spin combinations to obtain a whole set 
of bounds for various TMD PDFs. 



2.5.2.2 Wilson Line out to Infinity 

We can accommodate for a Wilson line of the form given in eq. (2.19) in the proof above 
by choosing a gauge G where the gauge fields vanish at infinity, e.g., the Feynman gauge. 
Then the transverse pieces are unity, and it is sufficient to define U[Ci t e] = U[£, oov + £] and 
M[C2.e] = M[oov, 0]. Obviously this choice satisfies the requirements above, with IA\C^^(\ = 
ZYfoou, 0], provided v has no functional dependence on t. Thus positivity of the respective 
quark-quark correlator is shown. 

The argument above has been presented without the inclusion of any subtraction factors that 
appear in the improved correlator definitions of section 2.4.4. It is clear, however, that a 
subtraction factor as part of the quark-quark correlator can destroy the positivity argument 
if it introduces an additional ^-dependence. 



2.5.2.3 Straight Wilson Line 

For the quark-quark correlator with a straight Wilson line U\C(\ = U[£, 0], it is not possible 
to find a prescription to split the Wilson lines in such a way that eq. (2.35) is satisfied, so we 
need to take a different route. Splitting the line "at infinity", we choose U[Ci^] = U[£, oo£], 
U[C2,e\ = U[ool, 0]. Let us introduce again quark fields with a "gauge history" 



(2.38) 
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The gauge invariance of the is guaranteed if we can assume that the gauge fields vanish 

at oo£ [RS65]. 4 This assumption is problematic for lightlike £. 5 Therefore, let us assume that 
£ is spacelike for the moment. The idea of gauge invariant fermion fields dates back to Dirac 
[Dir55]. Gauge invariant fermion fields with Wilson lines out to infinity have been discussed 
within QED, e.g., in Refs. [Man62, RS65], see also Ref. [JMY05] . 

We can now rewrite the first steps of our proof eq. (2.34) as 

i>[ rop ](^) = (n(p,s)\ *l/£) (r d )t r d ^,(o) \n(p,S)) = 

52 ( (n\ T d 9 qtt (£) \N(P, 5)) ) * ( (n\ T d 3^(0) \N(P, S)) 

n 

J2e i{p - pn > £ (n|T d <& q/ {0)\N(P,S)) 2 . (2.39) 

n 

Owing to the gauge invariant quark fields, we have been able to perform these steps without 
fixing a gauge up to now. However, in order to carry out the Fourier transform of $[ r ° P l(£) 
with respect to £, we need to show that the ^-dependence of the squared matrix element in 
the equation above is superficial. 

To this end, we now fix the gauge to the radial gauge G (also called Fock-Schwinger-gauge) 
[Foc37, Sch51, Cro80]. To obtain this gauge, we set W(£) =U[0,£] in eqns. (1.4)-(1.6). This 
procedure fixes the gauge completely up to a global color rotation, see also Refs. [Cro80, 
LW96]. The gauge field -A'(£) in the radial gauge satisfies = for any £, from which 

follows immediately that radial Wilson lines U'[£, 0] become unity for any £. Thus we obtain 

Y,e^ P - p ^ e ( G (n\T d q '(0)\N(P,S)) G )* . (2.40) 

n 

We would now like to continue the proof as in eq. (2.37). However, £ can become lightlike in 
the integral of the Fourier transform. If the integrand is regular for lightlike £, this does not 
cause any difficulties, because we can exclude null sets. If we make this assumption, eq. (2.37) 
is applicable, where now U[Cs] = 1. In this last step, translation invariance of the gauge 
condition G is not required. 

We have to acknowledge that the use of gauge invariant quark fields relies on non-trivial 
assumptions, which deserve to be checked again carefully in the future. 

2.5.3 Difficulties regarding the Interpretation of TMD PDFs 

Let us distinguish three different questions regarding the interpretation of TMD PDFs: 
1. Are we able to isolate features of the nucleon with TMD PDFs? 



4 To maintain this, valid gauge rotations W(x) of the gauge transformations eqns. (1.4)-(1.6) must become 
constant at foo. 

5 Keep in mind that for £ 2 = 0, the physical, Lorentz invariant distance (oo£) 2 = oo-0 is undefined. 



2.6 More TMD PDFs: Parametrization of $ [r ° P] 



31 



2. Are TMD PDFs universal, i.e., is the same set of distributions applicable for the de- 
scription of a multitude of scattering processes? 

3. Do TMD PDFs have a mathematical interpretation in terms of probability distributions? 

Concerns that we might have to answer the first question in the negative were raised when 
it became known that final (or initial) state interactions play an important role in parton 
distribution functions even at leading twist [BHM + 02] . It was argued that, due to the influence 
of the final state, parton distribution functions do not solely contain information encoded 
in the wave function of the nucleon. However, according to our present understanding, the 
relevant final state interactions are encoded in the gauge link, see, e.g. Ref. [Col02]. Therefore, 
it is not necessary to know the wave function of the final state in order to calculate parton 
distributions. Thus we may still think of TMD PDFs as properties specific to the nucleon. 

Regarding the question about universality, it has been argued that the same TMD PDFs 
describe both the Drell-Yan and the SIDIS process, apart from known sign changes due to 
the reflection of the gauge link [CM04]. The situation appears to be more complicated for 
other processes, see Ref. [BM08]. 

Let us now discuss the third question. The parton densities given by TMD PDFs must 
be positive and normalizable if we want to interpret them as probability distributions in a 
mathematical sense. In the previous section, we have discussed an approach to establish 
positivity of parton densities. Concerning normalizability, we encounter the same difficulties 
as already mentioned in section 2.3.2: In J dx f d 2 k± \x, k±; P, S), the fcj_-integral is 
undefined. The integral becomes finite if we restrict the integration range to < Q. 

Whether there is a more advantageous way to introduce a regularization scale is still a matter 
of ongoing research. With the Gaussian ansatz we will use in section 4.6, the fcj_-integral is 
finite without an explicit cutoff. Here the scale dependence is hidden in the ansatz, which we 
know to be inapplicable at large k±. 

In general, parton densities can only be interpreted as probability distributions within the 
context of an appropriate regularization and renormalization scheme and with respect to the 
corresponding scales. 

For all practical purposes in this work, we will give our results an interpretation in terms of 
quark densities intrinsic to the nucleon in the sense of probability distributions. The discussion 
above should make the reader aware of the more subtle issues regarding this point of view. 



2.6 More TMD PDFs: Parametrization of $I r ° P ] 

In the previous sections, fi(x,k±) served as an important example of a TMD PDF. Seven 
more of these profile functions are needed to describe spin dependent effects at leading twist. 
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2.6.1 Lorentz-Invariant Amplitudes 

In order to find out which TMD PDFs exist, the first step is a parametrization of the correlator 
(j>[r op ] i n terms of Lorentz-invariant amplitudes. Here we will briefly review this procedure, 
using a formulation that will be convenient for our lattice calculations and paying special 
attention to the role of the Wilson line. Let us look at a somewhat generalized quark-quark 
correlator that also allows spin-transitions of the nucleon: 

& r ° P \k,P,S,S') =/t04 e-* k - el 2 (N(P,S>)\ q(t)r*U[C t ] ?(0) \N(P,S)) . (2.41) 

We can parametrize this object in terms of Lorentz-invariant amplitudes that are independent 
of the nucleon spin vectors S and S'. This can be seen by rewriting the correlator in the form 

& r ° v \k,P,S,S') = l - U(P,S')Mro P (k,P;C e )U(P,S) . (2.42) 

Here U(P,S) are the Dirac spinors introduced in section A. 4 and A4r°p(k, P; Q) is a Dirac 
matrix. Notice that A4r°p{k, P; Cg) is not ^-dependent; rather, the appearance of Ci in the 
argument of M reminds us that there is a prescription how to connect the quark fields for any 
given quark separation I. The symmetry transformation properties of Air°p(k, P\ Ci) under 
Hermitian conjugation (f), parity (fp) and time reversal (T) are 

(t): [M T {k,P;C^ 
(*): M T (k,P;Ct) 
(T): [Mr(k,P;Ce)}* 

Here the bar indicates sign change of the spatial components of a four-vector, e.g., k = (k°, k) 
=4> k = (k°, —k). The path Cg is meant to be the reverse path of Q, i.e., Q(A) = Cg(l — A), 
and C is the charge conjugation matrix 7 2 7°. 

2.6.1.1 Straight Wilson Lines 

For the straight Wilson line U[Ce] = U[£,0], the Dirac matrix Mr°p(k, P;Ce) only depends 
on the four-vectors k and P. For a given r op , we now express A4r°p(k, P; Ci) as a linear 
combination of all Lorentz-covariant structures that can be formed from k and P, weighted 
with invariant amplitudes Ai(k 2 , k-P). We eliminate all structures that are incompatible with 
the transformation properties (f), (fp) and (T) given in eqns. (2.43)-(2.45). Note that the 
straight Wilson line fulfills U[Cg\ = U[C^i + £] = U[Cg) = U[—C_i\, so all three equations 
(2.43)-(2.45) provide useful constraints. Inserting the structures into eq. (2.42) and using the 
Gordon identities in eqns. (A.12)-(A.15), we find that some of the structures are redundant 
for the parametrization of & r ° P \k, P, S, S'), and we leave them out. We thus arrive at the 
parametrization given in eq. (B.l) in the appendix. In terms of <J?[ r ° P ](&:, P, S), the structures 



= 7° M lQ rt 7 o (k, P; C-t + t) 7° , (2.43) 
= 7°.M 7 or 7 o(£ ) P;Q) 7° , (2.44) 
= 7 5 C AW r . 7 s c (fc, P; -C~i) . (2.45) 
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thus obtained read 6 

<Z> [1] {k,P,S) = 2 TTiTv At , 
*M(fc, P, S) = 2 A 2 P M + 2 ,4 3 & M + 
2 



2 



m N 

+ 2 A 10 e^Saka + 



m N 



+ 2 A 9 e^S a P p 
2 A xl e^k a Pp(k ■ S) 



^r'](k,P,S) = -2m N A 6 S>" 
^(k,P,S) = [2i A 5 (k-S)} 



m N 



A 7 P^{k ■ S) 



m N 



A 8 k»(k ■ S) , 



(2.46) 



The structures above correspond to the parametrization worked out in Refs. [RS79, TM95, 
MT96]. Powers of have been inserted to make the A4 dimensionless. All amplitudes A± 
above are real valued because of the requirement A* = A{ following from the hermiticity con- 
straint (f). The structures proportional to A4, A§ and A\2 (highlighted by square brackets) 
change sign under time reversal (T), and are thus incompatible with eq. (2.45). An experi- 
mental measurement giving a non-zero value for such a T-odd amplitude can be interpreted 
as a clear sign that a non-trivial Wilson line is needed in the definition of the quark-quark 
correlator, as we will see in the next section. 



2.6.1.2 Wilson Lines out to Infinity 

Let us consider a Wilson line of the shape U[Ce] = U[Ci tV ] = U[£, r\v + £, rjv, 0] for rj — > 00 as in 
eq. (2.27) and, for 77 = 00, in eq. (2.20). This contour introduces an additional dependence on 
v, so the Lorentz-invariant amplitudes are now functions Ai(k 2 ,k-P,v-k,v 2 ,v-P). Only the 
directional information contained in v is relevant, so we may rescale v by and the 

amplitudes take the form 



A{ I k , k'P, 



v-k 



v-P 



v-k 



A i (k 2 ,k-P, r --,C 1 ,sgn(vP) 



\v-P\ ' \v-P\ 2 ' \v-P\J V" '" " ' \v-P\ 

Note that for v ~ n_, we have v-k j v-P ~ x and £ ~ 0. The Wilson line fulfills 
(t) : U[C- e , v +£}= U[£, r)v + £, V v, 0] = U[C t , v ] , 



(2.47) 



(fP): U[C lv ]=U[C tt€ ] 



(f): U[-C_ lv ]=U[C^\. 



(2.48) 



-v m 



from which we conclude that v transforms as (f): v — » v, (fp): v — * v, (T): v - 
eqns. (2.43)-(2.45). We recognize that time reversal converts a future-pointing Wilson line 
into a past-pointing Wilson line [Col02] . The arguments of the amplitudes Ai remain invariant 



6 Note that no additional amplitudes appear in the correlator $' r P \k, P, S, S') as compared to the "spin 
diagonal" correlator $f r ° P l (fc, P, S). 
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under application of (f) and (fp), but the third and the last argument change sign under (T), 
because v-k — ► — v-k = —v-k and v-P — > —v-P = —v-P. Obviously, time reversal (T) provides 
a relation between two subtypes of amplitudes, Ai{. . . , +1) and Ai(. . . , —1), rather than a 
relation that limits the number of allowed structures. This means that amplitudes such as 
the ones in square brackets in eq. (2.46) become possible [Col02]. 

A complete parametrization of <£ r ° P , which also takes the u-dependence of the Lorentz- 
covariant structures into account, has been published in Ref. [GMS05] and involves 32 ampli- 
tudes in total. In a lattice calculation with non-straight Wilson lines, the existence of all these 
amplitudes must be taken into account. Initial studies with a gauge link U[£, rjv + £, r]V, 0] on 
the lattice will be presented in section 5.2. 



2.6.2 From Amplitudes to TMD PDFs 

The next step is to rewrite the k~ integral eq. (2.6) for the different structures in (2.46) in 
Lorentz invariant form to isolate the profile functions. The procedure will be shown in more 
detail in sec. 4.1.2, when we address the Fourier transformed amplitudes. At leading twist, 
one obtains the following TMD PDFs [RS79, TM95, MT96, GMS05]: 



^ + \x,k ± ;P,S) = f 1 (x,k 2 J 



€± V k± * S± 3 f± ,V 1.2 

m N 



^ + ^(x,k ± ;P,S) = Ag 1L (x,k 2 ± ) + g 1T (x,k 2 j 



& a ^x,k ±] P,S) = e l [S 3 h 1T (x,k 2 ± ) + 



+ 



k l 

m N 



i 



m N 

f ij k- 
e ± K j 

m N 



(2.49) 
(2.50) 



Ahi L (x,k 2 ± ) + 



k± ■ S± 

m N 



hu<(x, fcj_) 



(2.51) 



where i,j = 1,2 are indices denoting transverse directions. The nucleon spin has been de- 
composed as in eq. (A. 8) in the appendix. Again, the structures in square brackets are T-odd 
and therefore not present for correlators implemented with a straight gauge link. Note that 
the integrated PDF corresponding to g\L is called g±, and that hi is the integrated PDF 
corresponding to h\x + {k\/2m? N )h^ T . The list of TMD PDFs at leading twist given above is 
complete, and does not change when the f-dependence is taken fully into account [GMS05]. 



2.7 Azimuthal Asymmetries from Transverse Momentum De- 
pendence 

Experimentally, we have access to the transverse momentum dependence of parton distri- 
butions through certain azimuthal asymmetries. Which asymmetries can be exploited to 
extract specific TMD PDFs is summarized in Ref. [BM98]. Below, we list some asymmetries 
of particular interest. 



2.7 Azimuthal Asymmetries from Transverse Momentum Dependence 



35 



2.7.1 The Cahn Effect 

Even in unpolarized measurements, we see evidence of the intrinsic transverse motion of quarks 
[Cah78]. The dependence of the SIDIS cross section on fa is governed by the fcj_-dependence 
of ^j^poiOcj k±; P, S) = fi(x,k^_). If the quarks carried no intrinsic transverse momentum, 

i.e., if ^J^po! (x, k±; P, S) = fi(x)5^ 2 \k±), then the dependence on fa would be trivial, c.f. 
eq. (43) in Ref. [Kot95]. 

Making the Ansatz 

1 / k 2 \ 

fl tq (x,k 2 ± ) = / 1|g (a?) -^ exp I J ( 2 - 52 ) 

motivated by the Fermi motion of partons in hadrons [Cah78, CT08], the authors of Ref. [A + 05] 
find that a value 

(fc3_)~ 1/2 = 0.5GeV (2.53) 

for the root mean square (RMS) transverse momentum agrees best with data from EMC 
[A+83, A+87b, A+91] and FNAL E665 [A+93a]. Based on the same ansatz, Ref. [CT08] 
addresses the question how changes under evolution of Q 2 . 

Lattice results related to the Cahn effect will be presented in section 4.6.5. 



2.7.2 T-Odd Effects 

The bracketed term in eq. (2.49) involves the so called Sivers function fyj<(x, kj_), introduced 
in Refs. [Siv90, Siv91]. This leading twist T-odd TMD PDF is proportional to k± x S± and 
thus describes the correlation of the intrinsic transverse momentum of unpolarized quarks 
with the transverse spin component of the nucleon. In SIDIS, the Sivers effect is accessible 
from the transverse spin asymmetry oc sin(^ — <j>s)- Therefore, a polarized nucleon target is 
required. Experimental information about the Sivers function has been obtained from Belle, 
HERMES and COMPASS data, see e.g., [D'A08]. Further experimental knowledge may also 
come from the PAX, PHENIX, RHIC and STAR experiments [C + 06b, C+06a]. 

Another leading twist T-odd TMD PDF, the Boer-Mulders function h±(x, k±), is accessible 
from unpolarized experiments by measuring a cos(2<^) asymmetry [BM98]. It can be inter- 
preted as a correlation proportional to k± x s± in the quark density, where s± represents the 
spin of a transversely polarized quark. Experimental results come from NA10 at CERN as 
well as E165 and E866/NuSea at FNAL, see, e.g., Ref. [ZLMS08, D'A08]. 

The existence of non-vanishing T-odd TMD PDFs, like the Sivers function, was ruled out 
at first [Col93]. Ideas that T-odd TMD PDFs may exist none the less due to final state 
interactions [BM98, BHS02a] could be cast into mathematical form, once the directional 
change of the Wilson line under time reversal had been discovered [Col02]. As mentioned 
before, the direction v of the gauge link switches from future-pointing to past-pointing when 
comparing SIDIS with the Drell-Yan process. Therefore the time reversal transformation (T) 
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transforms SIDIS TMD PDFs into Drell-Yan TMD PDFs. The T-even distributions are equal 
for the two processes, while the T-odd distributions change sign [Col02]. 

Obviously, T-odd TMD PDFs cannot be explored on the lattice with the straight Wilson line 
we use in chapter 4. With a staple-shaped Wilson line, we will get a first glimpse into the 
realm of T-odd contributions from the lattice in section 5.2.3. 



2.7.3 The Quark Density in the Polarized Nucleon from gir(x, k 2 ± ) 

In the pursuit of obtaining a picture of the quark density with respect to transverse momen- 
tum, let us write down the operator 



1 + A7 5 1 7+ + A7+7 5 

q\ +) W[Ci] ?(+)(o) = ^= mm] 7 2 7 7 • (2.54) 

This tells us that the density of longitudinally polarized quarks with helicity A is obtained 
from 

$ [(7 + +A7+7 5 )/2]( X)fc±;j p )<S ) = }-q>h + ]( Xjk± ;P,S) + ^ [,y+ ^ ] (x,k ± ;P,S) . (2.55) 
If the nucleon is polarized in transverse direction (A = 0), this becomes 



PTL(x,k ± ; S±, A) == \fi{x,k\) + ^ k± — g 1T (x,k 2 ± ) 
2 2 rriN 



l e ±ij k ±i S ±j ± 
2 rriN 



(2.56) 

where the bracketed term is T-odd (and vanishes for straight gauge links). We now recognize 
that g\x{ x i k\) parametrizes a correlation between quark helicity, transverse nucleon spin and 
transverse momentum, which introduces an axial asymmetry in the density. The structure 
with §it{x, k\) vanishes if we average ptl(%, k±; S±, A) over k±. Accordingly, an integrated 
PDF corresponding to q\t{x , k\) does not exist. Reference [BBHM00] gives the following 
bounds on gir{x, k\): 

\ (W) 2 + (/iV) 2 ) < (ft + 9il)(H ~ 9il) < Ui? , (2-57) 

from which follows 

k 



< l/il • (2-58) 
It is clear that this inequality must hold in order to guarantee positivity of eq. (2.56). 



giT — 

m N 



Experimentally, g\T(x,k\) is accessible from the azimuthal asymmetry cos((f)h — ^5) of the 
cross section. Its measurement requires longitudinally polarized leptons and a transversely 
polarized target [BM98]. 

In section (4.6.5), we will compute the first x-moment of g\T{x } k\) on the lattice using 
straight Wilson lines and will find that eq. (2.57) numerically also holds for the first Mellin 
moment. 



2.8 Model Predictions 



37 



2.8 Model Predictions 



TMD PDFs have also been addressed within models, see, e.g., the appendices of Ref. [MMG07]. 
To give an example, the authors of Ref. [JMR97] calculate the correlator eq. (2.5) within a 
spectator model (termed "scalar diquark model" in Ref. [MMG07]), and find 

f , ,2n AT2 {l-x) 2a - l {xm N + m) 2 + k\ 

h r(x,k I ) = JSIp, — s— , (2.59) 

Jl,HK ±J R lg7r 3 (^ + A |(x)) 2Q 

/ i2s (1 — x) 2 "" 1 xrriN + rn 

9lT , R (x,k ± ) = a R m N N R — 3 (fc a + > ( 2 - 60 ) 

where \r(x) = A 2 (l — x) + xM R — x{l — x)m 2 N and where the index R = a or R = s 
denotes the type of spectator. The TMD PDFs for individual quark flavors are obtained from 
fl,u = (3/2)/i )S + (l/2)/i )0 , fx,d = fi,a and analogously for g XT . At the given order, the 
antiquark TMD PDFs (ascribed to values x < 0) are zero. The authors of Ref. [JMR97] fix 
the constants Nr from the condition 

J dx J d 2 k ± h,R{x,k\) = 1 , (2.61) 

which ensures that the total number of up quarks is J dx f d 2 k± fx,u(x, k±) = 2 and the 
total number of down quarks is J dx j d 2 k± fx,d{ x i k±) = 1. Following this strategy and 
taking a typical set of values A = 0.5 GeV, m = 0, = 0.938 GeV, M s = 0.6 GeV, 
a s = 1, M a = 0.8 GeV, a a = -1/3 and a = 2 from Ref. [JMR97], we evaluate the first 
Mellin moment (see section 2.3.3) of the distributions numerically and plot them in Fig. 2.10. 
The rederived results in Ref. [MMG07] have a set to 1 and thus reproduce the perturbative 
tail fx{x,k\) ~ l/fcj_ at large fcj_ as mentioned in section 2.3.2. To allow for a qualitative 
comparison, we plot them in Fig. 2.10, leaving all parameters except for a unchanged. To be 
able to normalize them, we restrict the integration range to |fcj_| < 1 GeV in eq. (2.61). 
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Figure 2.10: Examples for predictions for the first Mellin moment from the scalar diquark 
model [JMR97] using the expressions and parameters specified in the text. The solid graphs 
were evaluated with a = 2 and are normalized to satisfy quark counting eq. (2.61). The 
dashed curves correspond to a = 1 and normalization with the cutoff 1 GeV. 



Chapter 3 

Nucleon Structure from Lattice 
QCD 

3.1 Basics of Lattice QCD 
3.1.1 The Path Integral formulation 

How do we calculate correlation functions within QCD? The path integral formulation of 
quantum field theory provides a very elegant prescription, and serves here as a starting point 
for the development of numerical methods. 

Suppose O is some expression in terms of quark fields q, q and gluon fields A, for example 
O := u(0) d(0) d(y) u(y) . According to the path integral formulation, we can calculate the 
vacuum expectation value (0| TO |0) of the corresponding time ordered operator TO from 

vv " JVA JVq JVq exp(tS QC D[g,g,A]) 

In the path integral above, q, q and A denote functions of space-time (the field configura- 
tions), and the expression O is a functional of them. The dynamics of the system arises from 
constructive interference in the integral over all field configurations; here T>q, T>q and DA 
represent appropriate integration measures. The action Sqcd is also a functional of q, q and 
A, and defined as 



/T /'oo 
dx° / d 3 x £ qCD [q,q,A)(x) . (3.2) 
-T J~oo 



The infinitesimal constant 9 > implements a tiny Wick rotation which makes sure that the 
state |0) in eq. (3.1) corresponds to the vacuum, see e.g. Ref. [PS95], chapter 9.2. 
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3.1.2 Path Integral in Euclidean Space 

For our purposes here, lattice quantum field theory serves as a tool to calculate matrix ele- 
ments in QCD numerically. One of the problems of a numerical evaluation of eq. (3.1) is the 
oscillatory term exp(i Sqcd[<Z, Q, A]). We therefore set 6 = it in eq. (3.2), or equivalently, we 
substitute 

xi = ix° = ixo , x\ = x 1 = —x\ , X2 = x 2 = —22 , x% = x 3 = —23 . (3-3) 

The new coordinates xp, are Euclidean: x^x^ = —xp,xp_ = — [x\x\ + 2322 + 2523 + 24X4). It 
is customary to define the Euclidean action Sq CD as 

*Sqcd[<7> Qi ^4.] = -iS Q cD[q,q,A]\ d=n 

= Jd 4 xl E q{x){$ + m q )q{x)+ l -F^ a {x)F-^ a {x)\ , (3.4) 

\q=u,d,s,... I 

where, in the last line, we have switched to Euclidean notation as explained in appendix A. 7. 
Equation (3.1) becomes 

, , = JVA JVq jVqO[q,q,A] exp(-Sg CD [g, q, A]) 

{ } " JVA jVq JVq exp(-^ CD [^,A]) 1 ' J 

The exponential term is now well-behaved, since <Sq CD [(/, g, ^4] is non-negative and real. Note 
that the 90° Wick rotation we have performed requires the analytic continuation of the fields 
q, q and A to Euclidean space. For the evaluation of 0[q, q, A], the fields q(x), q(x) and 
A{x) are in practice only available at Euclidean coordinates x. For example, the two-point 
correlation function obtained by setting O := u(0) d(0) d(x) u(x) can only be calculated for 
Euclidean separations x, since 24 is real. 



3.1.3 Discretization of Free Fermions 

By discretizing the action, we are able to replace the functional integrals in eq. (3.5) by a 
large number of ordinary integrals. Let us introduce a lattice {2(0), 2(1), £(2)> • • •} = L = aZ 4 
of points in four-dimensional Euclidean space with uniform lattice spacing a, and abbreviate 
vectors of length a along the Euclidean axes with fi = ae^. Consider the action of free quarks 
of a single flavor 

*.[ f „] = /A«,)(Vii + ™.)^). (3.0) 
A naively discretized version of eq. (3.6) is 

S%L%q\ =« 4 E IE fa) 7A ^ + A) ~ a q{X ~ A) + fa) m q q{x) ) , (3.7) 
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where dp has been replaced by a central difference. S\^ e depends on a countable number 
of field variables, namely q(xr \), q(x^), q(x^), q(xm), . . .. It is convenient to work with 
dimensionless variables on the lattice. Therefore, we will follow the convention that the quark 
fields on the lattice are rescaled according to q{x) — > a~ 3 ^ 2 q(x), and the masses are replaced 
by rh q = am q . The lattice action then reads 

i] = E f y, fa) 7m g(x+A) : 9(x " A) + fa) s fa) i . (■>■> 




Note that this action is bilinear in q and g, and can thus be written in the form 

SbL[q,q] = ^ x ) K (a,x),(f3,y) ipiy) . ( 3 - 9 ) 

where we have made Dirac indices a and (3 explicit, and where the Kr ax \tpy\ form a large 
matrix of coefficients. 

For given q and q, the discretized action S\^ c converges to the continuum action Sf ree in the 
limit a — > 0. Even though, this naive action Sj^ e does not provide the desired continuum 
limit in the path integral. Reasons for this will become clear in the following. 

3.1.4 The Fermion Doubling problem 

The dispersion relation that belongs to the lattice action Sffcc m ^ ne P rev i° us section has 
an unwanted feature: It has 15 additional energy minima at nonzero three-momentum. The 
problem is a consequence of the use of the central difference operation in eq. (3.8). The 
shortest wavelength that can be realized in one lattice direction is 2a. The central difference 
q{x + fi) — q{x — fi) vanishes for such a field configuration, because the central difference 
extends over 2 lattice spacings rather than just one. In other words, the central difference 
is blind to certain modes on the Brillouin zone. If we were to carry out a lattice calculation 
with the naive fermion action SfeL, all 16 low energy modes present in the action would be 
populated and would contribute to our measurements. We would have 16 fermion species 
instead of one. A number of techniques has been developed to deal with this problem. Here 
we list just three of them (For details, see, e.g., Ref. [Rot97, DD06].) 

• Wilson fermions: We can add a term 

4V^V^ * ^fa ~ A) - 2 fa) + fa + A) /o ln >, 
-a 2^2^arq{x) (3.10) 

xeL p,=i 

to the naive fermion action eq. (3.7). Here f is a dimensionless parameter. The expres- 
sion above is a discretized version of —{ar /2)q{x)dp j dpq{x). This term vanishes in the 
limit a — > 0. At finite lattice spacing it alters the dispersion relation, raising the energy 
of the spurious modes. Thus only the physical fermion species survives. The drawback 
of the Wilson action is that it does not exhibit chiral symmetry: When we set the quark 
masses m q to zero, QCD is invariant under the global transformation q — > exp(ie7g)g. 
At nonzero lattice spacing, the Wilson term destroys this invariance explicitly. 
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• Staggered fermions [KS75]: Consider again the naive fermion action eq. (3.8). It is 
invariant under translations of step size a along the lattice axes. It turns out that 
the fermion matrix Kt a ^ x \(p^ defined in eq. (3.9) can be block diagonalized into four 
blocks. This means that there are four independent groups of fermion species, without 
any mutual interactions between those groups. Therefore, we can simply remove three 
of those groups. The resulting action is invariant under translations of step size 2a, 
i.e., the unit cell of the action is of size (2a) 4 . This "staggered" action still exhibits 
16/4 = 4 fermion species, called tastes. Quark correlation functions calculated with the 
staggered action exhibit taste splittings, i.e., depending on the taste degrees of freedom 
that have been used to set up the correlator, the expectation value will be a little 
higher or lower. The taste splittings should vanish in the continuum limit a — > 0. The 
staggered action is invariant under a certain modified chiral transformation. A serious 
problem of staggered fermions is multiplicity: In fermionic loops, all four tastes give a 
contribution. To avoid this, one commonly makes use of the fourth root trick: In the 
Monte Carlo sampling step (see sections 3.1.8 and 3.1.9 below), one takes the fourth root 
of the fermion determinant. Whether it is guaranteed that this procedure reproduces 
QCD in the continuum limit is a matter of ongoing debate [Cre08, Gol08]. At present, 
we use staggered actions in spite of this open issue, because they are computationally 
exceptionally cheap, and thus permit us to explore parameter ranges and lattice sizes 
that would otherwise be inaccessible with present computing resources. 

• Domain wall fermions [Sha93]: This action (just as the overlap action) is able to establish 
a modified version of chiral symmetry on the lattice and is at the same time doubler 
free. The modified chiral symmetry transformation is q — ► exp[ie7g(l — (a/2fo)D)]q, 
where the operator D fulfills the Ginsparg- Wilson relation [GW82], namely {D,7g} = 
(a/fo)-D7g.D. Chiral symmetry on the lattice is desirable because it reduces operator 
mixing (see section 3.1.13 below), simplifies renormalization and is valuable for chiral 
extrapolation. The idea behind domain wall fermions is to separate left-handed and 
right-handed quarks spatially in an auxiliary fifth lattice dimension. The additional 
dimension comes at a considerable computational cost. The larger the lattice size L5 
in this auxiliary dimension is chosen, the more accurately chiral symmetry is fulfilled. 
Before correlation functions are evaluated, the lattice is projected onto the usual four 
dimensions. 

We note that the fermion actions above are still of the bilinear form of eq. (3.9). The calcula- 
tions in this work will be carried out using a hybrid action, with domain wall valence quarks 
on top of a staggered S6cl, clS explained in section 3.3 below. 

3.1.5 The Gauge Principle on the Lattice 

Upon discretization, some symmetries of the continuum theory are lost, but they are restored 
in the limit a —* 0. Concerning local gauge symmetry, it turns out that we can construct a 
discretized action that retains gauge invariance even at finite lattice spacing. Consider the free 
fermion action S^f in eq. (3.9). First, let us focus on quark bilinears involving neighboring 
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sites, such as q a [x) K(a,x),(0,x+ii) Q/3i x + A)- Under the gauge transformation eq. (1.4) and 
(1.5), this term becomes q a {x) W '(x) -^(a,x),G8,a:+A) W(x + p) qp(x + p). To make it gauge 
invariant, we introduce a new set of fields Un(x) of SU(3) color matrices, which transform 
according to 

U p (x) - %{x) = W{x) U^x) W\x + p) . (3.11) 

With appropriate insertions of the U fields, we can modify the quark bilinear terms in such a 
way that they become gauge invariant. For example, the terms 

q a (x) K {a)X)) ^ )X+(l) Up,{x) qp{x + jj) and q a (x + fi) K {a)X+(l) ^ )X) Uj^x) qp{x) 

are gauge invariant. The Un(x) are called link variables, and are depicted as lines connecting 
neighboring lattice sites, see Fig. 3.1. For later convenience, we introduce the notation 

Un(x) = U{x,x + p) , Ut(x) = U{x + p,x) . (3.12) 

Quark bilinears of non-neighboring sites require several insertions of adjacent link variables. 
For example, q a (x)Kr ajX \rp jX+ 2(i)U'( x ) x + P)U(x + fi, x + 2fi)qp(x + 2fi) is gauge invariant. 
Note that the prescription of rendering bilinears invariant is not unique; any connected path 
between the quark fields could be chosen. Finally, the resulting fermion action can be written 
in the form 

S l $%q,U] = ^2 Yl ?<*,i( x ) K (a,i,x),(p,j,y)[U] qp,j{y) = «LiJ K m\[U] <7|2J , (3-13) 

x,y£h a,/3 i,j |1J.|2| 

where we have made color indices i and j explicit. The products of link variables needed 
to maintain gauge invariance of the individual quark bilinears have been combined with our 
former coefficient matrix to form an object Kf a ^ x s (g^ y '\\(T\. On the right hand side, the 
indices and [2] each abbreviate a combination of a Dirac index, a color index and a lattice 
site. 

On the lattice, we get a rather intuitive understanding of the gauge principle. The quark field 
variable q a (x) is a color vector with respect to a frame of reference of color coordinates. Local 
gauge transformations allow us to rotate this frame of reference independently at each lattice 
site. The link variable Unix) tells us how the color frames of two neighboring lattice sites are 
rotated relative to each other. Thus they enable us to form expressions involving quark fields 
at different sites that remain unaffected by local color coordinate transformations, i.e., local 
gauge transformations. 

3.1.6 The Gauge Action on the Lattice 

Now we still need an action Sq [U] that determines dynamics of the link variables. We can 
give the link variables an interpretation in terms of the gluon fields An{x) of the continuum 
theory by writing 

U(x, x + fi)= exp (igaAp,(x)) . (3-14) 
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Making this identification, we request that S'(j t [J7] corresponds to the continuum gluonic part 
of the continuum action eq. (3.4) when formally taking the limit a — > 0: 

lim Sg[U] = [ d A x ]rF flPa (x)F fiDa (x) . (3.15) 

a^O J 4 

A simple action fullfilling this requirement is 

= £ EE^ 1 - u i»w) • ( 3 - 16 ) 

Here (3 = 6/g 2 is the coupling constant, tr c the trace over color indices, and the plaquette 
Ufiu{x) a rectangular loop of four link variables: 

Ufip(x) = U(x, x + ft) U(x + fx, x + p, + z>) U(x + p, + v,x + C>) U(x + v,x) . (3.17) 

Since Up,p(x) = Up^lx), the terms in the gauge action add up to a real, non- negative value. 1 

Together with the fermion action, we have now a lattice action S" lat [g, q, U] = S l f[q,q,U] + 
SgfiT]. Using this discrete action in the path integral eq. (3.5), we can replace the functional 
integrals by ordinary integrals: 

j VA J] (J dUfi(x) \ = JdU, Jvq^l[(J dq n } = j ' dq 

fj,= l. A 

and analogously for q. The path integral then reads 

/nS = JdU [dq JdqO[q,q,U] exp(-5 lat [g, g, ^) , , 

K ' JdUjdqJdqex V {-S^[q,q,U}) ' 1 ' 



3.1.7 Finite Volume 



We cannot simulate an infinite number of field variables. Therefore, we restrict ourselves to 
a finite lattice volume. We achieve this by introducing periodic (or antiperiodic) boundary 
conditions on the borders of a four-dimensional box. Typically, a lattice volume has three 
equal spatial side lengths of L lattice units and a temporal extent of T lattice units. A two- 
dimensional illustration of the periodic lattice is shown in Fig. 3.1. The number of integration 
variables in the path integral is now finite. 

In order to analyze nucleon structure, the nucleon should fit inside the lattice volume. From 
the point of view of chiral effective field theory, the nucleon is a core dressed by a cloud of pions. 
If the lattice box size is too small, our observables will be affected by pion mediated interactions 



As a side remark, consider the limit f3 — > oo. In this limit, all plaquettes are forced to unity. Consequently, 
any gauge link made up of link variables connecting two given lattice sites gives the same SU(3) c -matrix. In 
that sense, color space is "flat" , and we recover the free fermion action. 
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Figure 3.1: Illustration of the periodic lattice and its degrees of freedom. 
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of the nucleon with its mirror images on the periodically continued lattice. Therefore, the 
pion compton wavelength should be small compared to the lattice box size: < L. The 

lattices considered in this work feature m n L > 6, a size at which volume effects are usually 
smaller than statistical errors. In the following, we will make no effort to estimate systematic 
uncertainties from finite volume effects. 



3.1.8 Integrating out Fermions 

In the path integral formalism, the fermionic field variables qu and qu\ are Grassman variables, 
i.e., anticommuting numbers. The integrals over the Grassman degrees of freedom q and q in 
the path integral can be evaluated analytically. To this end, consider the generating functional 

Z[fj,r),U] = / dq j dq exp - ^ «|iJ -^[iJL 2 !^] q \A "XXiJ^iJ "XX^LiJ ■ ( 3 - 19 ) 

V L1J.L2J W LiJ / 

We complete the squares and perform the integrals over the Grassman variables to get 



Z[fj, v , U] = det (K[U\) exp [ fj^ Gyy [2J [U] m 

J1J-L2J 



(3.20) 



where the lattice quark propagator Cyy ^ [U] is the inverse of the fermion matrix: 

K [l] |2J [U] G [2} L3J [U] = <*[1J |3J • ( 3 - 21 ) 

The generating function allows us to integrate out the fermions for any expression 0[q, q, U] 
in the path integral: 



0[q,q, U]exp { - ^ q^ K^[U] q^ J = O 
L1J.L2J 



-,--,u 



ri=ri=0 



= det (K[U]) O 



d_ _d_ 
dr]' dfj 



exp [ ^77 L1J G [m [U]r) [2} 
LiJ 



Z[fj,v,U] 

= det (K[U]) 6[U] . 



»?=r?=0 



(3.22) 



With the help of the formula above, we are able to determine a representation 0[U] of O free 
of any fermionic variables. As a simple example, suppose we have 0[q, q, U] = 9|i|9|2J ■ Then, 
according to Grassmann algebra, we obtain 



Q[U] 



drj di] 



Z[rj,rj,U] 



r]=r)=0 



(3.23) 



Integrating out fermions amounts to forming all possible Wick contractions of quark fields 
in 0[q, q, U] and replacing each contracted pair by a propagator [U]. Finally, the path 

integral has a form suitable for numerical treatment: 



fdUO[U] det(K[U}) exp(-S>f[£/]) 
JdU det(K[U]) exp(-S^[U]) 



(3.24) 
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Note that up to now, we have restricted ourselves to a single quark flavor. For several 
quark flavors q = u,d,s, . . ., the fermion matrix K[U] = diag(ifW[C7], K^[U], . . .) is block 
diagonal, each block being a fermion matrix for a single flavor with components ^nj^j " 
Correspondingly, we can introduce lattice quark propagators GmL [U] for each flavor. 



3.1.9 Monte Carlo Calculations 



In eq. (3.24), we have to integrate over all independent link variables in our volume. For 
multidimensional integrals with very many dimensions, Monte Carlo techniques, such as the 
Metropolis algorithm, become very efficient. At each sampling step, the algorithm produces 
one gauge configuration U, an array of real numbers specifying the values of all the link 
variables on the lattice. In a simple Monte Carlo algorithm, each gauge configuration is a 
modification of its preceeding configuration. Not all configurations proposed by the random 
generator are accepted. The rejection policy is set up in such a way ("detailed balance") that 
the resulting sequence of gauge configurations Ejj, called the ensemble, contains an arbitrary 
configuration U with a probability proportional to the weight factor in the path integral : 

P(U e E v ) oc det(K[U]) exp(-Sjf *[C7]) . (3.25) 

Of course, a requirement for this importance sampling approach is that the weight factor 
det(if[Z7]) exp(— S^ft/]) is a real and positive number. This is why the Wick rotation to 
Euclidean space is essential. With an ensemble of N configurations at hand, we can now 
calculate correlation functions according to 

<°> 4^ d[u] ■ (3 - 26) 

ueEu 

For reasons of efficiency, modern lattice calculations typically combine Monte Carlo techniques 
with deterministic sampling algorithms, such as molecular dynamics. However, the basic 
definition of lattice expectation values eq. (3.26) remains unaffected from the particular choice 
of algorithm. 



3.1.10 Gauge Fixing 

The path integrals eqns. (3.2), (3.4), (3.18), (3.24) are all of the general form 



M.M^, (3,7, 



where subsumes all fields q, q, A and/or U and where the weight factor 9[4>] contains the 
exponentiated action (and, if q and q are integrated out, the fermion determinant). Gauge 
fixing can be achieved with a mapping f g which transforms a given field configuration 4> i n 
accordance with the gauge transformation rules eqns. (1.4)— (1.6), (3.11) into a new configu- 
ration (j)' which fulfills the gauge fixing condition g. To this end, the gauge rotation matrix 
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W(x) = W g [(j>](x) must be chosen appropriately for a given 0. The integration measure T><j) 
and the weight factor 9[<f>] are gauge invariant, in particular 0[<tf>] = 0[f g [(f>]]- Therefore, the 
gauge fixed version of the path integral reads 

, x - Mi m (3 28) 

= — jwem — ' ( } 

If O is gauge invariant, i.e., if 0[(p] = 0[f g [4>]], then gauge fixing has obviously no effect. If, 
however, O is not gauge invariant, then gauge fixing is mandatory to obtain a meaningful 
result. This will be the case in sections 4.4.5, 4.4.8 and 4.5, where we will use the Landau 
gauge fixing condition dp i Ap(x) = 0. This condition becomes ^2 X ^Re tr U(x,x + p) = min! 
on the lattice [DD06] . The Chroma executable has the corresponding minimization algorithm 
built in, and is able to convert a given gauge configuration U to Landau gauge, i.e., Chroma 
provides the mapping^. 2 



3.1.11 The lattice as a Regularization Scheme 

Consider an oscillation of the fields with an Euclidean wave vector k = (%, k%, fcj, kg). To 
the discretized action, wave numbers \kn\ > it /a are indistinguishable from oscillations with a 
corresponding wave number in the range \kp\ < ir/a (aliasing effect). To phrase it differently, 
any mode described by the degrees of freedom of the lattice path integral can be uniquely 
assigned to a momentum inside the first Brillouin zone, with —it /a < k^ < n/a. In that sense, 
the lattice imposes an ultraviolet momentum cutoff it /a. Moreover, since we are working with 
periodic boundary conditions in a box of size (aL) 3 x (aT), wave numbers are multiples of 
27r/aL and 2ir/aT, respectively. Thus the "resolution" of momenta is limited. As a whole, 
the lattice provides both an ultraviolet cuttoff (through the lattice spacing) and an infrared 
cutoff (through finite lattice volume). Consequently, all quantities are finite on the lattice. 

All quantities on the lattice are expressed in terms of dimensionless numbers 6(a). The corre- 
sponding quantity in physical units is obtained according to # lat (a) = a de 0(a), where dg is the 
length dimension of the quantity. Physical observables, e.g., hadron masses, are not sensitive 
to the behavior of the theory at very small length scales and converge to their continuum 
value: # lat (a) — > Q cont as a —> 0. There are other quantities which are sensitive the ultraviolet 
behavior of the theory. They are only well-defined with respect to a given renormalization 
scheme and renormalization scale or renormalization condition. The renormalization scheme 
provided by the lattice depends on the details of the lattice action used. It is therefore desir- 
able to translate # Iat (a) to a renormalization scheme in the continuum, such as MS. Important 
examples of renormalization scheme dependent quantities are the quark masses m q and the 
coupling constant g. 



2 The gauge rotations VK 9 [C/](a;) determined during the gauge fixing process can be used to produce gauge 
fixed lattice quark propagators later on. 
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3.1.12 Determining the Lattice Spacing and Setting Quark Masses 



Notice that after rescaling the fields, the lattice spacing a appears nowhere explicitly in the 
action. The lattice spacing is controlled by the coupling constants of the action, in particular 
by the lattice gauge coupling g, or rather f3, in eq. (3.16). To make contact to our physical 
world, we calculate some dimensionful observable on the lattice for a given value of (5 and 
compare to experimental results. This way, we can determine a. We may then adjust the 
value of until we are close to the desired lattice spacing. 

Since we are going to work with unphysical quark masses, the observable used to determine 
a should be largely independent of quark masses. A common choice is the static quark 
potential, which is then compared to phenomenological models describing the spectrum of 
heavy quarkonia, such as bb states. Some details of this procedure will become clear in section 
4.4.7, where we are going to calculate the static quark potential in order to renormalize the 
Wilson line. The lattice spacing determined this way is subject to statistical uncertainties of 
the measurement, and, most importantly, to systematic errors inherent to the method. 

The lattice quark masses m q also needs to be tuned. Modern lattice actions typically incor- 
porate the three lightest quarks u, d, and s as dynamical degrees of freedom. The strange 
quark mass m s can be set to an approximately physical value. The light quarks u and d 
are usually chosen degenerate, m u = = m u d. The most convenient observable to specify 
the light quark masses is the pion mass m n , because, unlike the "bare" quark masses m u d 
on the lattice, it needs no renormalization, yet is very sensitive to m n( j. In chiral effective 
field theory, the pion acquires mass only through the explicit breaking of chiral symmetry 
caused by nonzero quark masses, and it follows the Gell-Mann, Oakes and Renner relation 
[GMOR68] oc m U( i + C(m^ rf ). The proportionality constant in this relation compensates 
all renormalization scheme dependence of the quark masses. 

The computational effort in lattice simulations increases drastically for lower quark masses. 
Only very recently, advanced algorithms and machines made first attempts possible to go 
down to a realistic pion mass of around 140 MeV. However, the ensembles we are going to use 
for our exploratory calculations feature pion masses no less than 500 MeV. With input from 
such large pion masses, an extrapolation down to the physical pion mass is prone to exhibit 
a large unknown systematic error. It will be interesting to repeat our calculations at lower 
pion masses. 

Heavier quarks c, b and t do not appear in the lattice action, i.e., "dynamical effects" of heavy 
quarks are neglected. Due to the large energy required to create a heavy sea quark pair, such 
fluctuations can only exist for a very short time, smaller than the lattice spacing. Unless 
observables involving heavy valence quarks are considered, we can "integrate out" heavy 
quarks. Effectively, the existence of heavy quarks merely amounts to small adjustments of 
the coupling constants. 
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3.1.13 Operator Mixing 

For simplicity, let us ignore gauge fields for the moment. Consider the local operators 3 

O (a) = q(0) q(0) I ?(0 ) g(0 ) , (3.29) 

a 

Oiia) = m (?(-£) - 2g(0) + qm - 7 SCaHng . a £ q(0)d^d p q(0) + 0(a 3 ) . (3.30) 

' Taylor expand ^— - ' 

A* A 

The first operator Oq looks like the mass operator in the action. Switching to dimensionful 
quark fields and dividing by the cell volume a 4 ( "rescaling" ) , we recognize the familiar X/a 
behavior of the mass term. The second operator looks like the Wilson term eq. (3.10) and 
the "naive continuum limit" obtained from a Taylor expansion of the fields reveals a relation 
to the Laplace operator in the continuum. The reason why the Taylor expansion is naive 
is that there is no one-to-one correspondence between the field configurations on the lattice 
and continuous fields. This is different from the classical picture, where the lattice can be 
chosen fine enough to give a rather accurate description of the smooth continuous field. It 
turns out the Laplace operator appearing on the right hand side of eq. (3.30) is only part of 
the continuum operator that corresponds to 2 . When we increase the lattice spacing the 
operator 02(a) on the fine lattice will have a representation in terms of O 2 (o!) and Oo(a') 
(and other operators) on the coarse lattice: 

02(a) = Z 20 (a, a') Oo(o') + Z 22 (a, a') 2 (a') + ... . (3.31) 

The Zij are matching factors. Operator mixing happens quite generally when switching the 
renormalization scheme or changing the renormalization scale. The continuum representation 
of 2 can therefore be specified in the form 

2 ( a ) - Z20(a '^ g(0) q(0) +ai 22 (a,rf^ 9 -(0)^#) + . . . (3.32) 



and the matching factors now depend on the lattice action, the lattice spacing a, the 
renormalization scheme in the continuum and the corresponding renormalization scale fi. 
The mass-like contribution to the continuum representation of 02(a) is potentially enhanced 
by a factor a~ 2 with respect to the second derivative term. 4 

An operator can mix with any other operator that has the same symmetry transformation 
properties under the symmetries of the action. An important symmetry of the Euclidean 
continuum action that restricts the number of operators that can mix is the 0(4) rotational 
invariance. The lattice is not rotationally invariant, but there is a remnant of the 0(4) 



3 In the context of path integrals, they are not operators. We will call them operators none the less - for 
convenience, and because of their obvious correspondence to operators in the Hamiltonian formalism. 

4 Indeed, for the Wilson fermion action, this means that the Wilson term eq. (3.10) contributes significantly 
to the quark mass through an additive renormalization. 
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symmetry, the hypercubic group H(4): The action is invariant under permutation of the axes 
and under reflections: 



Here (fxi, /J-3, ^4) is a permutation of (1,2,3,4). Since H(4) is less restrictive than 0(4), 
mixing patterns on the lattice are often much more complicated than in the continuum, see, 
e.g., Ref. [G+96, G+05a]. 

3.1.14 Action Improvement 

We have a lot of freedom in setting up the discretized action as an approximation to the 
continuum action. For example, loops of link variables other than the plaquette also yield 
the field strength tensor in the continuum limit. The systematic method to design optimized 
actions is called Symanzik improvement [Sym83a, Sym83b]. In essence, the design principle 
works as follows: All local lattice operators invariant under the desired symmetries of our 
action can possibly be part of the action. From these, one selects an appropriate finite 
subset of operators 0\ at . The improved action is a linear combination of these operators, 
with coefficients chosen in such a way that spurious operators in the corresponding "effective" 
continuum action cancel up to a certain order in a. However, as we have seen in the previous 
section, the coefficients appearing in the continuum action cannot be determined from a Taylor 
expansion of the lattice operators. Rather, lattice perturbation theory or non-perturbative 
methods must be used. For more information on action improvement, see, e.g., Ref. [DD06]. 

3.1.15 Link Smearing: HYP Blocking 

The gauge configurations determined from the Monte Carlo sampling contain a lot of noise, 
in the form of short distance, high momentum fluctuations. Since we are usually interested in 
the long distance behavior of the theory, it can be useful to suppress this noise. This can be 
done by one or several smearing steps. Smeared gauge links, also termed fat links, can be used 
during the production of gauge configurations as an ingredient to an improved fermion action, 
or can be used in subsequent steps of the analysis. The smearing technology applied in this 
work is a single step of "hypercubic blocking" , or HYP smearing [HK01] , which is implemented 
ready to use in the Chroma library. The fundamental operation of HYP smearing is an APE 
smearing step [A + 87a], where each smeared link variable is formed from the unsmeared link 
variable and an admixture of staple-shaped gauge links: 



U A pe{x, x + jx) = Proj su(3) ^(1 - o)U{x,x + p) 

+ f Yl U(x,x + u)U(x + u,x + P + fi)U(x + u + p,,x + p,)) . (3.34) 



The projection to SU(3) ensures that the new link variables are again unitary. For HYP block- 
ing, three iterations of APE-like steps are carried out in such a way that the final fat link 
variable receives contributions from links no farther away than the edges of the surrounding 



(x 1; x%, £3, x 4 ) > (rbx^j , ix^ 2 , ±2^3 , ix^ 4 ) . 



(3.33) 
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hypercubes. This way, the smearing operation remains local. Long range properties (univer- 
sality class of the theory, spectrum, etc.) remain unaffected. A disadvantage is the influence 
on short distance behavior, e.g., of the static quark potential [DD06], compare section 4.4.7 
and Refs. [A + 87a, HK01]. The smearing strength of HYP blocking is controlled by three 
parameters. The ones used by LHPC and throughout this work are a.\ = 0.75, «2 = 0.6 and 
a 3 = 0.3. 



3.2 Nucleon Matrix Elements 



This section will describe the method used to calculate properties of the nucleon on the 
lattice. Since we have used the Chroma programming library [EJ05] and numerical input from 
the Lattice Hadron Physics Collaboration (LHPC), we have taken over their conventions and 
techniques. 



3.2.1 Baryon Sources and Sinks 

We need to create nucleons on the lattice, although we do not know the precise wave function 
of the nucleon. The first step is to define a suitable "interpolating field" , a gauge invariant 
combination of quark and gluon fields with the same quantum numbers as the nucleon. We 
shall use 

1 



B a (t,P) 



'L 3 



E< 



-iPx 



e a b c Uaa(x,t) (u{ (x,t) T diq d c (x,t) 



(3.35) 



as a nucleon sink on the lattice. Here we have made the Dirac index a and the color indices 
a,b,c explicit. The nucleon source is located on a time slice t in Euclidean space and has 
three-momentum P. Rather than making the obvious choice r diq := C75 = 747275, LHPC 
uses r diq := (775(1 + 74). The non-relativistic projection can reduce the number of Dirac 
indices that need to be processed. The adjoint expression B a (t, P) serves as a nucleon source. 
The quark fields q(x,t) appearing in eq. (3.35) are replaced by "extended quark fields", a 
gauge covariant superposition of quark fields located in the vicinity of (x,t). This technique 
of source and sink smearing increases the overlap of the nucleon interpolating field with the 
nucleon wave function. For details, see Refs. [Gus90, AJG + 91, A + 93b, D + 02]. For later 
convenience, we write source and sink as 



B a (t,P) 



B, 



m E 



L 3 



iPx -ndiq _ ~j — 

e e L1J L2J L3J 1 |3J |2J U HJ d L2J U L 3 J 



*.|1J.|2J,|3| 



t (t,p) = ^= £ 



L 3 



— iPx -ndiq j 

6 £ L4JL5JL6J 1 |5J |6J M L4 U ® d Y 



cc,L4j,|5j,|6j 



|1J,|2J,|3J at (x,t) 
Dirac |lj =a 

L4J,|5),|6J at (x,t) 
Dirac |4|=c* 



(3.36) 



where it is understood that e|iJ|2J|3J ac * s on ^ m c °l° r space and rj^j on ^y m *he space of 
Dirac indices. 
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3.2.2 The Nucleon Two-Point Function 



Two-point correlation functions on the lattice are used for the calculation of hadron masses 
and will serve us as normalization factors in the computation of nucleon matrix elements. 
The nucleon two-point correlation function is introduced as 



/3a 



(3.37) 



where the Dirac matrix T 2pt = (1 + 73)(1 +i"/^^)/2 used by LHPC is again a non-relativistic 
projection, combined with a spin projection. The latter is not needed in principle, but is 
statistically advantageous in combination with three-point functions. Inserting the source 
and sink, we get 



C 2pt {t,P) 



»,LiJ,..., 



-iPy 



-pdiq p2pt ^ 



; 14 L5J L6J 1 L5j L6J 1 [ij L4I ftl I 2 ! |3J 1 L3j L2| 



x («L4I U L5J d ® U IH d M U L3|) 



L1J,L2|,L3) at (aJsrcW) 
|4J,|5J,|6J at (y,t+t src ) 



(3.38) 



where we have already exploited the translation invariance of the expectation value. In the 
second line of the above expression, we form Wick contractions in order to integrate out the 
fermions: 



^ L 4jM|ij lt|_5jU[_3j <%jd|2| 



u L4R3j "KiJ d^d^} 



LlJ,L2|,L3| at (a; src ,t src ) 
L4J,L5J,L6J at (y,t+tsvc) 



(3.39) 



3.2.3 Point-To-All Propagators 

Each of the contracted pairs in eq. (3.39) becomes a propagator, see section 3.1.8. Since 
we have rh u = rhd, u- and d-quark propagators are identical. Moreover, notice that all 
propagators are attached to the source location at x src = (x SIC ,t STC ). Therefore, it is sufficient 
to prepare a single type of point-to-all propagator from the inversion 

K m\ M G m M = L3J I L3J at , arc ■ (3.40) 

The restriction to a single source location is crucial for the computational feasibility. This way, 
we only need to perform ~ N 2 x 4 2 inversions of type (sparse matrix) x (vector) = (vector) for 
each lattice configuration. The non-relativistic projection reduces the number of inversions 
further. 

The inversion above produces forward propagators Gpj [3| = ?|2|?|3j > where the second index [3j 
is fixed to the source location. We can use the relation (based on the so called "75-hermiticity" ) 

= G L3jL 2j = (75^75)12] L3J ( 3 - 41 ) 

to obtain a backward propagator, which has the first argument fixed at the source location. 
This relationship works for Wilson-type fermions, including Domain Wall fermions, and can 
save the large computational costs of further inversions. 
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Once the propagators are prepared, the nucleon two-point function can be evaluated in Chroma 
very efficiently, as illustrated in the following code snippet 5 

LatticePropagator diquark = 

quarkContractl3( prop * BaryonSpinMats : : Cg5NR() , 

BaryonSpinMats: :Cg5NR() * prop ); 
return trace( BaryonSpinMats :: TmixedO * traceColor( prop * diquark ) ) + 
traceColor( traceSpin( BaryonSpinMats :: TmixedO * prop ) * 
traceSpin( diquark ) ) ; 

After projecting out the desired momentum _P, this yields C 2pt (i, P) for all t. 
3.2.4 Transfer Matrix Formalism 

In order to be able to interpret the correlation functions, we switch to the transfer matrix 
formalism [Liis77]. In section 3.1.1 we claimed that a lattice expectation value is related 
to the vacuum expectation value of O. On a periodic lattice of finite extent T in Euclidean 
time, this is not precisely true. Rather, we obtain a trace. Specifically, for the nucleon 
two-point function, we obtain if < t < T: 

(B a (t,P)B (Q,P)) = fT "* B a (0,P) ¥ %(0,P) \n) 

n 

= TV T T - 1 B a (0, P) f 1 2^(0, P) . (3.42) 

Here {\n)} is a complete set of states which are normalized to unity. The transfer matrix T 
can be formally written in terms of the Hamilton operator H. If the states |n) are eigenstates 
of the Hamiltonian with energy eigenvalue E n , we may write 

f * = exp {~Ht\ = \ n ) e~ Ent (n| (3.43) 

n 

and substitute the expression on the right at each occurrence of T. Suppose now that T—t^>t 
is large, such that we can substitute f T - t « |0) in eq. (3.42). If t is also large, such 
that excited states are suppressed, the propagating state between the nucleon interpolating 
operators must be the nucleon \N(P, S)). Lighter hadrons or the vacuum are not allowed due 
to the quantum numbers of the interpolating operators. So we get 

— ^ P -Ept ^_ 

(B a (t,P)Bp(0,P))*J2(n\B a (0,P)\N(P,S)) -— (N(P,S)\B (p,P) \Q) . (3.44) 

s p 

The nucleon energy in the denominator appears due to the Lorentz covariant normalization 
of the nucleon state, (N(P, S)\N(P', S')) = 2Ep 5ppi 5sS'- Introducing overlaps 

Bp(0, P) \N(P, S)} = ^Z N {P) U a (P, S) , 

(N(P,S)\%(0,P)\n) = ^/Z N (P)*U a (P,S) (3.45) 



5 Similar code can be found in Chroma: : Baryon2PtContractions : :sigma2pt. 
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the nucleon two-point function eq. (3.37) finally reads 

C*P\t, P) = \ z ^ p )\ e -E P t tro r 2 P t y u{p ^ s) ^ (P) s) 
2Ep is 

1 y n 2E P lhpc 1 y ,l Ep 

(3.46) 

To be precise, we should mention that there is another state that can contribute: the antipar- 
ticle of the nucleon's parity partner. This state comes with the overlaps 

(tt\Bp(Q,P) \N'(P,S)) = y / Z N ,(P)V a (P,S) , 

(N'(P,S)\%(0,P)\Q) = ^Z N ,(P)*V a (P,S) (3.47) 

and thus yields a contribution proportional to trD[F 2pt J2 S V(P, S) V(P, S)] in the two-point 
function. Through our specific choice of r 2pt , this contribution will be largely suppressed, in 
particular for small momenta P. We only see the effect of the nucleon parity partner in the 
two-point correlator when we approach i src from the "wrong" side, t < t src , see Fig. 3.4. 

Now, from the slope of In C 2pt (t, P) = —Ept + const, we can immediately read off the energy 
Ep of the nucleon state. For the special case P = 0, we obtain the nucleon mass m^- 

3.2.5 The Nucleon Three-Point Function 

We want to be able to calculate nucleon matrix elements as they appear in eq. (2.41) on the 
lattice [MS89, WDL92]. The "operators" we encounter are of the form 

Or° P , q (C) = q(l[C])T°P U[C] q(0) , (3.48) 

where £[C] = C{1) — C(0). The quark type q can be u or d. We use a lattice field combination 
of the form 



O^P i9 (C lat ; z) = q{£[C l&t ] + z) r op Z/ at [C lat + z] 

^op ;/lat[y-ilat 



to create analogous operators on the lattice. We have introduced an offset z, which leaves 
the matrix element in eq. (2.5) invariant. Here U lat [C lat + z] is a product of connected link 
variables starting at position z and ending at position £[C lat ] + z, following the convention 
eq. (A. 25). Contracting with the color indices of the two quark fields, it renders Opo P gauge 
invariant. The choice of the link path C lat will be discussed later. However, let us demand 
from the start that the link path remains on a single time slice, or in other words, that it does 
not contain any link variables in 4-direction. We now consider the three-point function 

C^ 9 (T,P;C lat ) = i EE 1 ?? P ) 0^ Pi(? (C lat ;z,r) Bp(t SIC ,P)) . (3.50) 

z /3a 
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B B B B 




Figure 3.2: Illustration of diagrams contributing to the nucleon three-point function on the 
lattice, here for an operator probing ci-quarks. Each tick line connecting quark fields represents 
a lattice propagator. The lower two diagrams are neglected in our calculations. 



Let us calculate the three-point function for (i-quarks in the operator. Substituting source, 
sink and operator, and exploiting translation invariance, we get 



C^(r,P;C lat ) = £ 
w,LiJ,..., 



e y e L 4JL5JI 



-pdiq p3pt ( 



[5] [6] 1 [1JI4J WJL 3 ! 1 |3| L 2 J 



r L7j L8J ^ lat t c ' at + (*» r )l L7J L8J 



LlJ , L2J , L3J at x src 

L4] , L5J , L6J at (y,t ank ) 

jV] at (z,t), LSJ at (£[C lat ]+z,r) 



(3.51) 



There are four possible contractions, as depicted in Fig. 3.2. Two of the contractions are 
disconnected, i.e., the quark propagators connect the operator with itself rather than with 
the nucleon. The calculation of disconnected contributions is computationally challenging due 
to a bad signal to noise ratio, and very expensive, because it requires all-to-all propagators. 
Meeting these challenges is a hot topic of lattice QCD today. However, throughout this work, 
we shall ignore disconnected diagrams, assuming that their effect is small. Note that discon- 
nected contributions cancel in isovector quantities, i.e., when we take the linear combination 
of operators 0^ p u _ d = O 1 ^ u - 0$ p d . 



3.2 Nucleoli Matrix Elements 



57 




Figure 3.3: Calculation scheme for nucleon three-point functions on the lattice, here for an 
operator probing d-quarks. Three lattice quark propagators of each of the two contributing di- 
agrams in Fig. 3.2 are combined into a single sequential propagator, displayed as the dark blob. 
To evaluate the three-point function, we calculate the product of the point-to-point propa- 
gator between the <i-quarks, the sequential propagator (for a fixed nucleon three-momentum 
-P), and a product of link variables from the underlying configuration U lat . 



3.2.6 Sequential Propagators 

Already for the connected diagrams, it looks as though we need an all-to-all propagator: The 

quark propagator cZ[y] has neither of its indices [6] and [7J attached to a fixed lattice site. 
However, in this particular case, there is a trick to avoid all-to-all propagators [WDL92], see 
also the appendix of Ref. [D + 02]. The two connected contributions can be brought into the 
form 

c r°p,d( T > p > clat ) = 73 r °7J |8J ^ lat ^ lat + T )l 17J |8J 

L2J,L6J,L7J,L8J L z 



X 



[2| at #src 

|7| at (z,t), |8J at (£[C lat ]+z,r) ' 

(3.52) 



Fd(x BTC , P, £ sn k) [2J [ 



LiJ , L3J , L4J , L5J 



-iP y , 



-pdiq p3pt ^ 



^diq 



L5JL6J 1 [1J|4J e LIJ L2] L3J 1 |3J L2J 



// 1 1 - 1 

X « L4J U |5J u Yi\ u \z\ 



r 



U L4 j U]^ U[ij U[_ 3 j / 



|1J,|2J,|3J at x SIC 
L4J,L5J,L6J at (y,i snk ) 



(3.53) 
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We now introduce the sequential propagator 

Sd(x SIC ,P, W)|2||7| = y^-^C^srCi-P^spk)^!^ ^|6|^|7| | [2| at x BIC • (3.54) 

|6J 

Multiplying the fermion matrix ^[71 [gj from the right, we obtain 

Sd(x STC , P, t snk ) |2j |7j ^^{gj = F d (x src , P, t sn k) |2J |9J | L2J at x BTC • (3.55) 

L7J 

Obviously, the sequential propagator can be computed from a matrix-vector inversion, sim- 
ilar as the point-to-all propagator is calculated in eq. (3.40), just with the sequential source 
Fd(x STC , P, t sn k) on the right hand side instead of a Kronecker-Delta. Sequential propagators 
for operators O^op u with u quarks can be set up in an analogous way, only the expression for 
the sequential source eq. (3.53) differs. The three-point function now simply reads 

C^Jt,P;C^) = £ ^T°V m U^[C^ + (z,T)] im 
L2J , m , L8J * 

X (S q (x SIC , P, t snk ) [2 j L7 j d [s] d [2] » ® * [8J at (£[c lat ]+z>T) • (3.56) 

Notice that the sequential propagator appears like a backward propagator in the expression 
above, with the first argument fixed at the source location. It is common practice to store 
sequential propagators in the format of forward propagators, so we must apply eq. (3.41) 
before we can use them for our purposes: 

int G5 = Ns*Ns-l; 

LatticePropagator GF = Gamma (opGamma) * fwprop; 
LatticePropagator seqpropbackw = Gamma(G5) *seqprop*Gamma(G5) ; 
return locallnnerProduct ( seqpropbackw, GF ); 

The Chroma code above 6 calculates the three-point function (up to the sum over z) for a gauge 
path of zero length, using a forward propagator fwprop and a sequential propagator seqprop 
loaded from file in advance. The integer value opGamma allows us to select any Dirac matrix 
r°P listed in Table B.l. 

The sequential source technique described above enables us to calculate nucleon matrix ele- 
ments for any operator insertion 0^ P u / d - Note, however, that separate sequential propagators 
need to be calculated for each of the quark flavors u and d, for each source-sink separation in 
Euclidean time, for each spin projection r 3pt , and, most importantly, for each nucleon sink 
momentum P. 



3.2.7 Assembling Gauge Paths 

The path C lat appearing in eq. (3.49) is a sequence of adjacent lattice sites starting at the 
origin. As a compact notation, we specify link paths by a sequence of "moves" - integers in 



Comparable code can be found in the file lib/meas/hadron/BuildingBlocks_w . cc of Chroma. 
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the range {—3 ... — 1,1... 3}. A positive number [i represents a shift ft, a negative number 
— n describes a shift —fi in the opposite direction. The paths that our program trans3pt 
processes can be specified in an XML control file. For the assembly of link paths, we used 
the algorithm in lib/meas/hadron/BuildingBlocks_w. cc. In essence, the propagator fwprop in 
the code snippet above is replaced by a propagator shft_fwprop, which is initialized with a 
forward propagator and then iteratively accumulates link paths and shift operations through 

int mu = spacedir (linkdir , j_decay) ; // "mu = llinkdirl" 
if (linkdir > 0) { // forward shift 

tmpprop = adj ( u[ mu ] ) * shft_f wprop; 

shft_fwprop = shift ( tmpprop, BACKWARD, mu ); 

} 

else { // backward shift 

tmpprop = shift ( shft_f wprop, FORWARD, mu ); 

shft_f wprop = u[ mu ] * tmpprop; 

} 

The program trans3pt automatically saves some computer time if a link path is an extension 
of the previous link pattern. 



3.2.8 Matrix Elements from Ratios 

After application of the transfer matrix formalism to the three-point function, the dominant 
contribution for <C r — t STC *C £ S nk — ^src "C T is 

C^ q (r,P;C^)^ \Z N (P)\ e -^(W-W) ^ U(P, S) T** U(P, S>) 

S,S' ^ p ' 

x (N(P,S')\ l-J2O^ P ^;(z,0)) \N(P,S)) , (3.57) 

where we have again used the fact that the contribution from the parity partner of the 
nucleon is suppressed for our choice of T 3pt . The sum over the offset z increases statistics in 
our computation. The unknown overlaps and the exponential can be cancelled by forming a 
suitable ratio with the two-point function: 



R F op tq {P;C lat ) = plateau 



C 2pt (t sn k — i src , P) 



£ 9F T'jr£tT U ^ (P ' g/) < °ffp, g (C lat ;0) MP,S)) . (3.58) 

^ 2E P tr D (-iP + m N ) j 

Here plateau [•] indicates that we are, in principle, referring to the limit r— t src — > oo, t S nk — r ~~ > 
oo, T — £ sn k + t src — * oo. In practice, we have to live with an approximation, of course. Plotting 
the ratio for a sufficiently large, fixed source-sink separation as a function of r, we expect the 
formation of a plateau between source and sink (Example plots will be shown in Fig. 3.5). 
We extract the value of Rr°p, q (P; C lat ) from this plateau region, within which we may assume 
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that contributions from excitations of the nucleon have decayed to a negligible level. A larger 
source-sink separation t src — i sn k reduces contaminations from excitations, but worsens the 
signal to noise ratio. 

At this point, it is necessary to establish a relation to renormalized continuum operators O ren . 
In general, the continuum operator of interest can be expressed as a linear combination of 
lattice operators Op<J p „(C lat ). The ratio Ro^(P) for such a renormalized operator O ren is 
then itself a linear combination of ratios Rr°p, q (P; C lat ), and provides access to the desired 
continuum matrix elements: 

^ 2E P tr D (-iP + m N ) ) 

We can now proceed analogous to section 2.6.1, writing the matrix element in the form 

(N{P, S')\ Icn \N{P, S)) = U{P, S') M ™{P) U{P, S) , (3.60) 

where A4o reu (P) is a Dirac matrix. Inserting this into eq. (3.58), and performing the sum over 
spins, we finally obtain 

p tr D {.M rcn(P) (-iJf> + m N ) T 3 ^ (-j# + m N )} 

Ro ren \P) = r n i , m ^ • (3.61) 

With the continuum parametrization M.o iou (P) a t hand, we gain access to the desired spin 
channels. We will make use of this master formula in section 4.1.2 and in equation (5.25). 
Sample ratio plots will be shown in Fig. 3.5. 



3.3 The Lattices of the MILC and LHPC Collaborations 

The computationally most intensive steps of a lattice simulation are the generation of config- 
urations and the calculation of propagators. In this work, we have used configurations and 
propagators calculated by the MILC and LHPC collaboration. For the purposes of these ex- 
ploratory studies, we have not selected the largest and most realistic lattices available today. 
Rather, we have chosen lattices of moderate size with pion masses of around ~ 600 MeV, 
which have allowed us to perform the analysis with reasonable statistics at a computational 
expense still manageable on the small PC cluster of our theory department. 

The gauge configurations we have used are listed in Table 3.1. They have been generated 
by the MILC collaboration [B + 01, A + 04, B + 07]. They feature 2+1 dynamical quarks, with 
the strange quark mass fixed to an approximate physical value. The gauge action is one-loop 
Symanzik improved (see, e.g., [ADL + 95, Bis05] and references in [B + 01]), and the AsqTad 
fermion action is of the staggered type, built using fat links and with order a 2 discretization 
errors removed at tree level, see Refs. [OTS99, Lep99] and references in [B + 01]. Note that 
the ratio rh U) d/rh s is constant for the last four lattices listed in Table 3.1. These ensembles 
serve us to study renormalization of the Wilson line operator in section 4.4. 
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ensemble-ID 


a (fm) 




rh s 


10/5 


?3 v T 
Li XI 


m DWF 


coarse-m050 


0.1181(17)(27) 


0.05 


0.05 


6.85 


20 3 x 64 


797(02) (29) 


coarse-m030 


0.1190(16)(27) 


0.03 


0.05 


6.81 


20 3 x 64 


621 (02) (22) 


coarse-m020 


0.1196(15)(27) 


0.02 


0.05 


6.79 


20 3 x 64 


514(02)(18) 


fine 


0.0854(12)(19) 


0.0124 


0.031 


7.11 


28 3 x 96 




superfine 


0.0601(08)(14)* 


0.0072 


0.018 


7.48 


48 3 x 144 




supercoarse 


0.1765(32)(38)* 


0.0328 


0.082 


6.485 


16 3 x 48 





Table 3.1: Lattice parameters of the MILC gauge configurations [B + 01, A + 04, B + 07] used 
in this work. The lattice spacing a in physical units is determined using r±, as explained 
in the text. The first error quoted for a includes statistical errors and the fit uncertainties 
specified for n, the second error stems from the systematic uncertainties of r\. The values 
marked with an asterisk * are preliminary [TD08]. The last column lists the pion masses as 
determined with the LHPC propagators with domain wall valence fermions [H + 08]. The first 
error is statistical, the second error comes from the conversion to physical units using a as 
quoted in the table. Note that the pion masses quoted here in physical units differ slightly 
from those listed in Ref. [H + 08], because we use a different strategy to fix the lattice spacing, 
see footnote 7. 



The numbers for the lattice spacing a in Table 3.1 come from the "smoothed" values r\ja as in 
Ref. [A + 04]. As conversion factor to physical units we always use the continuum extrapolated 
value r*i = 0.317(7) sys (3)fit fm given in Ref. [A + 04] for all ensembles 7 . 

We have used the first three lattices of Table 3.1 to calculate nucleon matrix elements, with 
propagators and sequential propagators from the LHPC collaboration. These propagators 
have been produced within a "hybrid action" approach: The calculation of nucleon matrix 
elements with staggered fermions is feasible but turns out to be rather complicated due to the 
unphysical taste degrees of freedom. Instead, the approach of the LHPC collaboration is to 
carry out the inversions eqns. (3.40) and (3.55) in five dimensions with a domain wall fermion 
matrix ifP|^[j7]. Even so, the gauge background U needed as input comes from the staggered 
simulations done by the MILC collaboration. The resulting propagators describe "domain 
wall valence quarks on a staggered sea", i.e. effects of virtual quark loops are included in the 
staggered formalism. The hybrid action approach is a suitable compromise that enables us to 
have chirally invariant, doubler-free valence fermions on a computationally affordable gauge 
background. For a physically meaningfull setup, the valence quark masses should be identical 
to the quark masses in the staggered sea. To accomplish this, LHPC tunes the valence quark 
mass parameters until their pion mass agrees with the mass of the lightest state in the pion 
taste multiplet ("Goldstone pion") directly obtained with the MILC staggered action [H + 08]. 



7 In contrast, the LHPC collaboration typically uses values for a obtained without the continuum ex- 
trapolation of ri but rather from a direct comparison of lattice results to the splittings of the T spectrum 
[WDG + 04, A + 04], which gives a = 0.124 fm on the coarse and a = 0.087 fm on the fine lattices. However, at 
present, the lattice spacings for the superfine and supercoarse lattices are only available in terms of r%. 
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The pion masses after this adjustment are displayed in the last column of Table 3.1. 

To reduce computational costs further, the lattices are chopped into two halves of temporal 
extent T/2 = 32. Only every sixth configuration, and alternating temporal halves are taken, 
reducing autocorrelations to an undetectable level. Noise is reduced by application of HYP- 
smearing to the gauge configurations before the inversions are performed, see section 3.1.15. 
The sequential propagators produced by LHPC are available for sink momenta P = and 
P = (—1,0,0) x 2n/L. The latter corresponds to \P\ ~ 500 MeV and is the lowest non-trivial 
momentum on these lattices. The source-sink separation is fixed to i S nk — ^src = 10. 

We should remark here that the point-to- all propagators generated by LHPC are by default 
calculated with smearing at the source. The sequential propagators have a smeared nucleon 
source and a smeared nucleon sink; only the free index is point-like. The two-point function 
in the ratio Rr°p, q in eq. (3.58) must be compatible with this setup, i.e., source and sink must 
be smeared. Therefore we first need to prepare smeared-smeared point-to-all propagators. 
This functionality is provided by the Chroma executable. 



3.4 Techniques of Statistical Error Estimation 

Powerful theoretical frameworks are available to deal with statistical errors in Monte Carlo 
data [Efr82, DD06]. For brevity, we restrict ourselves to a rather informal motivation of the 
prescriptions we have applied. 

The evaluation of n correlation functions d\ [U] , . . . , O n [U] for each of the N gauge configu- 
rations U in our ensemble yields samples {i = l..n, k = 1..N) from which we compute 
our observables. For each i, the samples sn~ scatter around some "true" expectation value Sj. 
Any of the final observables 9 we want to compute can be expressed as a function of the Sj. 
The "true" value is thus 9 = 0(si, . . . , s n ). Of course, the Si are not available to us. However, 
we can calculate an approximate value 9 = 0(s\, . . . , s n ) from the sample means 



1 N 

*i = j7 ^ Sik ■ (3.62) 



k=l 



We also need an estimate of the uncertainty of 9. Calculating the standard error of the sample 
mean from the scattering of 9 according to 



0-0 



^ N(N 



1 N 

T ^ ) jy<e(s lk ,...,s nk )-9) 2 (3.63) 



would be problematic, because the experience large fluctuations, much larger than the 
expected fluctuations of the Sj around the Sj. Since 9 can behave highly non-linear for such 
large fluctuations, we would introduce a large uncontrolled bias. The idea of resampling 
methods is to construct new samples Sjm that experience fluctuations of similar magnitude 
and distribution as the s,- around the h. 
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In general, we use the jackknife sampling method by default. In some cases, we opt for the 
bootstrap method, in particular when input samples from different ensembles (with different 
lengths N) need to be studied in a global analysis. 



3.4.1 Jackknife Sampling 

Initial jackknife samples are obtained from sample means with one input value missing. Those 
are used to form jackknife samples of 9: 

1 N 

s i(k) = N _i E Si i ' ^(fc) = ®( s i(k),- ■ ■ , s n (k)) • (3.64) 

3=1 

First of all, jackknife sampling [Que56, Tuk58] is a method to reduce bias in the final estimate 
9 . Let us interpret our input samples as random variables S^. The fluctuations in the mean 
values Si scale with iV -1 . Thus, in terms of expectation values E[-], it seems appropriate to 
expand the bias of 9 = 9(S%, . . . , S n ) in powers of A r_1 . The same thing can be done using 
jackknife samples instead of sample means: 

E[0(S!, ...,S n )}= 9(h, ...J n ) + ^ + (Jp) , (3.65) 

E[9(S 1M , S„ i(fc) )] = 9(h, . . . , Sn) + + O ( , N ^ 1)3 ) V k ( 3 - 66 ) 

with some coefficient a. This analysis immediately shows that the bias corrected jackknife 
estimate 

0Jack _ g _ J- ( 9 {k) - 9 ) (3.67) 

k 

has remaining bias of order N~ 2 . Moreover, the variances of the jackknife samples s^f.) are a 
factor (N — l) -2 smaller than those of the sik- This tells us how to modify eq. (3.63) in order 
to obtain a viable estimate of the standard deviation of 9, i.e., the jackknife error: 



^.Jack 



\ 



N-l N 



N 

k=l 




3.4.2 Bootstrap Sampling 

The bootstrap method [Efr79, Efr82] allows us to form an ensemble of samples whose size 
M is independent of the number of configurations N. In this work, we choose M = 100 or 
M = 1000. Larger values of M give more accurate error estimates but can be computationally 
costly. The first step of the resampling procedure is to draw N x M random integers a k i in 
the range 1..N. The bootstrap samples are then formed according to 

1 N 

Si ( l ) = at E Siakl ' =0( s l(i)>'--> J VO) ' (3.69) 

fc=i 
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In the sum above, some samples Sij may appear multiple times, others are left out. In order 
to handle correlations correctly, it is important that one and the same set of indices ctki is 
used for each i. The distribution of the s^m is an approximation of the distribution of the 
Sj. Likewise, we expect the bootstrap samples 9n\ to be distributed in a similar way as 9 is 
distributed. To determine an error interval, we simply sort the 9m according to their size, 
resulting in a sequence On^ < 0o 2 ) < . . . < 9n M y Let us throw away 15.8% of the entries at 
the start of this sequence, and 15.8% at the end of the sequence. The remaining sequence has 
a length of about 68.3% of M, and the lower and upper values define an error interval with 
confidence level of about 68.3% ("lcr"). 

Bootstrap sampling is particularly useful when combining data of several ensembles e, each 
with a different number of configurations N e . Then each channel i of our input data Sik 
has been obtained from a certain ensemble e^. Data coming from different ensembles are 
uncorrelated, so we may determine a new set of indices a e kl for each ensemble. Finally, we 
obtain resampled data of uniform length M just as in eq. (3.69), using the ofiS as random 
indices. 

3.4.3 Overdetermined Systems, Fits and Correlations 

In practice, the extraction of an observable 9 often involves solving an overdetermined system 
of equations. Typically, these equations are of the form Jj = fi(a\, . . . , a m ), with m < 
n parameters dj. The target observable 9 can be expressed as a function 9(a±, . . . , a m ). 
Obviously, given approximate input values s%, . . . , s n , we cannot find parameters that satisfy 
all n equations at once. Throughout this work, we will stick to the following strategy: 

The parameters dj are determined by minimizing numerically 

n 

X 2 = ^2 w i (/i(ai, • • • ,a m ) - Si) 2 . (3.70) 
i=i 

Here the choice of the weights Wi involves some arbitrariness. If the Sj are largely uncorrelated, 
a sensible choice is Wi = l/(Asj) 2 , where Asj is the jackknife/bootstrap error of Si itself. If 
the correlations among the are strong, it would be appropriate to modify eq. (3.70) such 
that the correlation matrix can be taken into account. Since the estimated correlation matrix 
is inaccurate and often close to singular, this approach would require special precautions. For 
this exploratory study, we will not follow this strategy, but rather adjust the Wi by hand 
according to simple, problem specific criteria. As a consequence, the minimized value of x 2 
does not have the strict statistical interpretation as it would have in a proper x 2 minimization. 

What is important to us is that the minimization yields values for the parameters at, from 
which we can determine 9 = 9{a\, . . . , a m ). Regardless of the difficulties with the choice of 
the Wi, this prescription implements 9 as a consistent estimator (see, e.g., [EDJ + 71]). 

We repeat the minimization eq. (3.70) with the entire set of jackknife/bootstrap resampled 
values Sj(n as input. This can be sped up if we pass the a±, . . . , a m as initial guess to the min- 
imization routine. From the jackknife/bootstrap estimates 9^ thus obtained, we determine 
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error intervals as described in the sections above. These uncertainties are appropriate for the 
estimator we have used. 

Because of correlations among the input data, the final estimate determined from our fit is 
statistically not optimal. However, the error bars we calculate with the resampling method 
take this into account. Their reliability is not compromised by our negligence of correlations 
in the fit procedure. 

3.4.4 Autocorrelations 

In the Monte Carlo simulation, the gauge configurations are determined from a Markov pro- 
cess, i.e., new configurations are random modifications of the preceding ones. Therefore, 
we expect correlations of a sample sn~ with preceeding samples Si k—i, Si,fc-2> • • •• Common 
techniques to estimate autocorrelations are binning or blocking, or the determination of an 
integrated autocorrelation time from jackknife samples 8 , see, e.g., [Jan02, DD06]. The LHPC 
collaboration has taken precautions to minimize autocorrelations (see section 3.3) and has not 
found autocorrelations of a detectable level in their observables. Therefore, for this study, we 
assume that autocorrelations can be ignored. 

3.5 Simulation Setup and Extraction Procedure 

Simulations of nucleon observables are carried out on the coarse lattices, where LHPC propa- 
gators are available. Remember that LHPC applies HYP smearing to the configurations and 
chops them in two halves of temporal extent T = 32. The point-to- all propagators have their 
fixed index (their "source") at t STC = 10. These propagators are the only ingredient needed to 
calculate the nucleon two-point function, which we display in Fig. 3.4. 

The LHPC sequential propagators we use have the nucleon sink at t s i n ^ = 20. The typical 
appearance of the three-point function is shown in Fig. 3.5. The oscillations of the correlator 
close to the source are an artifact that can appear when using domain wall fermions. To apply 
the formalism of section 3.2.8 we must evaluate the ratio Cpop (r, P; C lat )/ C 2pt (t sn k — t svc , P) 
at time slices r far enough away from source and sink (and the chopping boundary). We 
assume here that this is the case for the time slices f = 14, 15, 16. We obtain the final result 
for the ratio Rr°p(P; C lat ) from an average over the values from these three time slices. 9 
Table 3.2 lists the number of configurations we have used on the different ensembles. 



8 Reliable estimates of autocorrelation times require a large number N of configurations, ideally well beyond 
1000. 

9 There are more sophisticated procedures available to extract the ratio value. These take the effect of excited 
states (and possibly the oscillating states) into account, see, e.g., Ref. [R + 07]. Using a more sophisticated 
approach, one may obtain slightly different ratio values and more realistic error bars. Since the primary goal 
of this work is not a high precision analysis, we stick to the simple averaging procedure. 
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number of configurations 


ensemble-ID 


m T MeV 


HYP-smeared link 


unsmeared link 


coarse-m050 


pa 800 


425 


135 


coarse-m030 


ps 600 


563 




coarse-m020 


ps 500 


478 





Table 3.2: Number of configurations used for our calculations of nucleon three- and two-point 
functions. 
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Figure 3.4: Two-point function for the nucleon on the coarse-m050 ensemble. The circled data 
points have been used to fit exponential decays. The fit weights have been chosen according to 
the statistical errors of the individual data points. The decay to the right of the source belongs 
to the nucleon state. The nucleon mass thus extracted is tun = 1.649(08) s t a t(59) a GeV. The 
data points close to the chopping boundary should not be taken into account. 
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Figure 3.5: Plateaus plots for the real part of the ratio Cf^_ rf (Y, P, C lat )/C 2pt (t snk - t STC , P) 

on the coarse-m050 lattice for the link path C lat given by the moves 1,2,1,1,2,1,1,2,1 as 
depicted in Fig. 4.1. The gauge link was evaluated on the HYP-smeared lattices. The quark 
separation is \£\ ~ 6.7a ~ 0.8 fm. We plot the ratio for the two available nucleon momenta: 
(a) P = 0, (b) P = (-1,0,0) x 2vr/L. The plateau value R TiyU _ d (P; C lat ) is extracted from 
the three circled points and is displayed as a horizontal error band. 
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Chapter 4 

TMD PDFs with Straight Links 



In this chapter we present first results for TMD PDFs calculated on the lattice, albeit with 
a straight gauge link between the gauge fields. As discussed in section 2.4, this choice of 
the Wilson line does not correspond to the situation in scattering experiments. Thus the 
TMD PDFs we obtain cannot serve as quantitative predictions for TMD PDFs used in the 
description of SIDIS or Drell-Yan experiments. Nevertheless, a straight Wilson line appears 
to be a good starting point, for several reasons: 

• The simplest, most obvious way to render the bilocal operator gauge invariant is the 
straight Wilson line. From a purely theoretical point of view, such operators are ap- 
pealing in themselves. 

• We can study a discretized version of the straight gauge link and its properties directly 
on the lattice. 

• The parametrization of the matrix elements in terms of Lorentz invariant amplitudes is 
considerably simpler, see sections 2.6.1. 

• To our present understanding, TMD PDFs with straight links have a probability inter- 
pretation, see section 2.5.2.3. 

The fundamental object of interest in this chapter is the quark-quark correlator with a straight 
Wilson line 1 

i>[ r °%, P, S) = \ (N(P, 5)1 q(£) T°p U[£, 0] g (0) \N(P, S)) 
2 „ / 

= O r o P (£) 

= l - U(P, S') A4ro P (^, P) U(P, S) . (4.1) 
We shall work directly with a discretized version of Or°p {(■) ■ 



Here and in the following, we omit the index for the quark type q. 
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4.1 Parametrization 



4.1.1 Parametrization of the ^-dependent Correlator 

We have already discussed a parametrization of the quark-quark correlator in section 2.6.1. 
There the correlator and amplitudes were fc-dependent. For our lattice calculation, we need 
an analogous parametrization of the ^-dependent quantities. The symmetry transformation 
properties of the matrices Mr°p(£, P) read 



(t): 

(*): 
(T): 



Mt(£,P)\ =7° M j0rtl o(-£,P) 7° , 

M r (e,P)=7° M 7 o r ^(l,P) 7° , 
M T (£,P)]* =j 5 C M ch5ri5 c(-lP) C f ' 



(4.2) 
(4.3) 
(4.4) 



It turns out that we can find the structures compliant with the symmetry constraints by 
replacing k with im 2 N l in eq. (2.46) 2 . We thus obtain 



$ m (£,P,S) 
^ ] (£,P,S) 

& a ^(£,P,S) 
^° ] (£,P,S) 



2rriN A\ , 

2 A 2 P^ + 2im N 2 A 3 t + hi m N A 12 S v P a 
2im N A 4 {Pn v - P v £")] + 2 Ag e^ al3 S a P p 
+ 2im N 2 Aio e^Saip -2m 2 N A n e^£ a Pp(£ ■ S) , 
-2 m N A 6 S 11 - 2i m N A 7 P»{£ ■ S) + 2m N 3 A s ■ S) , 
-2m N 2 A 5 (£-S) 



(4.5) 



Again, the structures with A4, A§ and A12, highlighted with square brackets, are T-odd 
and vanish for our analysis with straight links. The amplitudes A{(£, £-P) appearing in the 
parametrization above are not constrained to be real valued like their momentum dependent 
counterparts Ai{k,k-P). Instead, hermiticity (f) now leads to the following relation: 



M£ 2 ,£-P) =A l (£ 2 ,-£-P). 



(4.6) 



4.1.2 Amplitudes Ai from the Lattice 

We now have a continuum parametrization of the ^-dependent matrix elements at hand, and 
can evaluate our master ratio formula eq. (3.61): 



pren / r> 



tr D {7Wro P (£,P) {-iff> + m N ) T 3pt {-i$> + m N )} 
2E P trD^P^-i^ + mjv)} 



(4.7) 



2 and likewise in the M structures eq. (B.l). For completness, the expressions for the Mr°p are given in 
eq. (B.2) in the appendix. 
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For J\Ar°p(£, P) we take the expressions listed in eq. (B.2), and for r 2pt and r 3pt we use the 
LHPC conventions. For any Dirac matrix r op , the equation above enables us to express the 
ratio -Rp p(-P,£) as a combination of amplitudes A{. We give a list of these expressions for 
16 basic Dirac matrices r op in Table B.l in the appendix. The relations in this table form 
a system of 16 equations, which we can solve for the individual amplitudes Ai. However, for 
this exploratory study, we have only picked a few simple channels, for which we do not need 
to disentangle the amplitudes: 

. 2A 2 = , (4.8) 
• 2A 6 = i£"L for gauge links with £ 3 = , (4.9) 

171 jy 13,5 ' ' 

For the value E(P)/mN needed to extract Aq, we use the continuum dispersion relation 



E(P) ^Jrn% + (2ir/L)2 

— = . (4.11) 

UlN TUN 

Here rh^ is determined as described in section 3.2.2 and illustrated in Fig. 3.4. We obtain for 
the ratio above 1.04942(47) sta t on the coarse-m050 lattice, and 1.0713(12) stat on the coarse- 
m020 lattice. Of course, as an alternative, we can also extract E{P) directly using a two-point 
function. The value thus obtained for E{P) /m{N) has a larger statistical error, and agrees 
with the value extracted using eq. (4.11) almost within error bars. 

Note that we will prefer to show results for mjsrA-j rather than the dimensionless amplitude A-j. 
The calculation of the latter would involve an explicit factor mjv, which introduces additional 
systematic errors and would complicate a chiral extrapolation. For the same reason, we will 
present results for the corresponding TMD PDF gyj* in terms of gir/mN or giT\k±\/iriN. 



4.1.3 From Amplitudes to TMD PDFs 



How are the amplitudes Aj,(£, P) related to TMD PDFs? Combining eq. (2.5) and eq. (2.6), 
we find that we need to perform the Fourier transformation 



& r ° P \x,k ± ;P,S) 



dl_ d 2 £ ± ei{ -e-k+ 



£+=0 



(4.12) 



2vr (2vr) 2 

where k + = xP + , as usual. We rewrite this integral in terms of the Lorentz-invariant quan- 
tities £ 2 and i-P. For £ + = 0, we have 



Thus we get 

^ {x , k± . P ,S) 



■P = tP 1 



■P) 



-£± ■ £± < . 



-i(£-P)x 



2tt 



£+=0 



(4.13) 
(4.14) 
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We substitute £± ■ k± = \£± \ \k±\ cos(0) and use 



/poo . . r2n -i rco pI'K 

d 2 £± = J d (7^2) ^/^£ 2 J d8=-J d(-£ 2 ) J d6 



Finally, we get 







p+ 


I'' 




dO 


Jo 


2^ 



-i{l-P)x 



2tt 



00 d{-£ 2 ) 
n 2(2vr) 



Jo 2?r 



£+=0 



(4.15) 



(4.16) 



Since ^ r ° P \i, P, S) is composed of amplitudes Ai, which are functions of £ 2 and £ ■ P only, we 
shall always be able to perform the ^-integral analytically. In the most general case, this can 
be done as shown in appendix B.1.3. Thus we are left with integrals over the Lorentz-invariant 
quantities £ 2 and £ • P only. 



4.1.3.1 Unpolarized case: f\{x,k\) 

In the straight link case, fi(x, fc^_) = ^ + \x, k±; P, S) (compare eq (2.49)). With £ + = 0, we 
get from eq. (4.5) ^ + \£, P, S)\ l+=0 = 2A 2 P + , which we substitute into eq. (4.16): 



h(x,kl) =^(x,k r ,P,S) 



d{- 



2tt 



cie 



i_ r d(£-p) i[e . P)x 

P+ J 2vr J 2(2vr) J Q 2vr 

2vr J 2(2tt) J q 2tt 



-e 2 \h±\ cos 9 2 j± 2 p^ 



d(£-p) e - iie .p )x r d{-e_ 

2tt J 2(2tt) 



e- |fe ± i cose 
JoCv^lfeil) 2i 2 • 



(4.17) 



Here Jo is a Bessel function of the first kind. With the notation introduced in appendix B.1.3, 
we simply write 



^(x,k r ,P,S) = T 0fi (x,\k ± \) [2A 2 
4.1.3.2 Axial vector operator: giz,(x, k±) and giT(x,k±) 



(4.18) 



Let us analyze ^^^(x, k±). The structure we read off from eq. (4.5) for the contribution 
from I 6 to $[t + ^ 5 ] is: 

- 2m N A G S + = -2 A A G P + , (4.19) 

where we have used the spin conventions eq. (A. 8) in the appendix. Then the corresponding 
contribution to 7 1 is 



^ 5] (x, k ± ; P, S) = A T ,o(x, \k ± \) 



-2A 6 



(4.20) 
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For A 7 , we have 



/ P + 

-2im N A 7 P + U-S)\ e+ _ n = (-2im N A 7 P+) At £ ± ■ S 

\ raw 



Applying the formalism of appendix B.1.3, we obtain 



^ + / ] (x,k ± ;P,S)=AT lfi (x,\k^ 



-2iA 



k±-S ± 

m N 



Tq,i(x, \k±\ 



•7 



2iAn 



Amplitude A& drops out since i + = 0. Comparing to eq. (2.50), we get 



g 1L (x,k 2 ± ) = T ,oO, |fci| 
gi T (x,k 2 ± ) = T i(x, I fei| 



-2A 6 
2i A 7 



+ Ti }0 (x, \k_ 



-2iA-, 



(4.21) 



(4.22) 



(4.23) 



4.1.3.3 Tensor operator 

We now analyze <&[ <7 * \x, fei; P, S), for a transverse index i = 1,2. Inserting £ + = and 
.Pi = into the corresponding line in eq. (4.5), remembering that A^ vanishes for straight 



links, and using e l+a P 



g a+ eP + g /3+ e' l f, we get 



&^l(£, P,S)=2 A 9 P + e^S a - 2i m N A 10 A P + eflp 

P+ 
m N 

We now apply the formalism of appendix B.1.3 and compare with eq. (2.51). Then 



2m 2 N A u P + e i ?£ a [ At i\ S } 



hyr{x, |fci|) = T 0i0 (x, |fei| 
hi L (x, |fei|) = T 0jl (x, \k±\ 
hi T (x, |fci|) = To >2 (x, |fci| 



2A 9 
-2i Aio 
2 in 



+ Ti,i(x,|fci|) 



-2 A 



ii 



(4.24) 



(4.25) 



4.2 The Discretized Non-Local Operator 

To realize the non-local operator Or°p(£) on the lattice, we approximate the Wilson line 
between the quark fields by a product of connected link variables along a lattice path C], at = 
(I, a;( ra_1 ), x( n_2 ) . . . , x^ 1 ), 0). We end up with a lattice operator just as in eq. (3.49): 

O&JC^O) = g(£)T op W lat [C] at ] q(0) = q(£)U(£,X^) ■■■ U(x^,0)q(0) . (4.26) 



If the Wilson line does not coincide with one of the lattice axes, i.e., if it is at an oblique angle, 
we approximate the straight line by a step-like path as illustrated in Fig. 4.1. To this end 
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Figure 4.1: Example of a step-like link path: The straight gauge link in the continuum with 
£ = (6, 3, 0) (dashed line) is represented as a product of link variables Un in the directions 
fi = 1,2, 1,1, 2, 1,1, 2,1. 



we programmed a Mathematica function resembling the Bresenham algorithm [Brc65], which 
lays out the path through lattice sites with minimal distance to the straight line connection. 

We interpret £ as a finite, physical separation of the quark fields. In the continuum limit a — > 0, 
the number of link variables required to form the gauge link goes to infinity: I = tja — > oo. 
The fine-grained details of the discretization prescription C\ at for the link path should not 
matter, as long as the discretized link can be made to lie arbitrarily close to the straight line 
by choosing a small enough. If this is maintained, the continuum limit of our lattice gauge 
link becomes a unitary path ordered product of gauge fields located arbitrarily close to a 
straight line. In that sense the continuum limit of our lattice operator Opo P q (C lat ;0) is the 
continuum operator Oy°p{£)- 

Obviously there are many ways we can realize a lattice gauge path close to a straight line. First 
of all, this raises the question whether the continuum limit is really unique and well-defined. It 
turns out that the continuum limit of our operator involves a power divergence, proportional 
to the lattice scale 1/a. The divergence and the respective ambiguity of the continuum limit 
must be adequately renormalized, see section 4.4. Together with the renormalization of the 
quark fields, we will end up with a relation of the form 



O r T %{£) = Z-\-£ 2 ; a) O^C^; 0) . (4.27) 



We will usually assemble our gauge links on a HYP smeared ensemble, so the individual links 
are in fact fat links. Effectively, the link we construct is a unitarized linear combination of 
many paths. The HYP smeared gauge links turn out to exhibit almost perfect rotational 
invariance with respect to £, indicating that they closely resemble a continuum operator. 
Without smearing, discretization errors appear more pronounced. 



4.3 First Observations on the Lattice 
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4.3 First Observations on the Lattice 

Let us focus on the channel r op = 73, and ignore the renormalization factor Z(£ 2 ; a) of 
eq. (4.27) for the moment. From eq. (4.8) we get from the unrenormalized ratio 

2Af 1Ten (l 2 ,£-P) := R T4 (P, 4 at ) . (4.28) 

As mentioned before, the amplitude A2 is related to the unpolarized TMD PDF fi(x, kj_). 
Let us begin our explorations on the coarse-m050 lattice, where we get the best statistics. 
Looking at the sample plateau plots shown in Fig. 3.5, we see a clean signal and surprisingly 
good statistics considering the rather long gauge link of about 0.8 fm. 



4.3.1 The Nucleon at Rest on the Lattice 

If we use the sequential propagators with the nucleon at rest, P = 0, we can only access the 
amplitude A% nren (£ 2 , l-P) at £-P = 0. Then the results we obtain from the ratio R ni {P,Cf i ) 
should only depend on the distance £ 2 = —£ 2 between the quarks; the residual dependence on 
the details of the link discretization C\ at should be small. 

To check this, we generate an initial set of link paths with the algorithm mentionend in the 
previous section. Firstly, we take link paths aligned with the lattice axes, with a length up to 
\£\ = 20 (We have studied link lengths up to \£\ = 40. Links on the axes that are longer than 
our lattice size L = 20 overlap with their periodic mirror images. We find that the signal 
remains compatible with for such very long links). Also, we cover quark separation vectors 
£ in the first quadrant of the (£\, ^-plane, as well as two quadrants in the (£\, £3)-plane for 
quark separations up to \£\ < 8. Moreover, we pick certain longer links in the first octant in 
order to cover more ^-values up to \£\ < 15. Last but not least, we pick some oblique links 
with integer lengths. The choice is illustrated in Fig. 4.2a. 

Calculating the ratio R Ti {P, C] at ) with these gauge paths on the coarse-m050 ensemble, we 
obtain Fig. 4.3. For the link paths evaluated on the HYP smeared lattice, we find an almost 
perfectly smooth dependence on £ 2 , even though we include step-like gauge paths. Obviously, 
the operator constructed from fat links exhibits rotational symmetry to a good approximation. 
This corroborates the notion that the discretized operator with an extended gauge link is an 
approximation to the straight link continuum operator. The situation is a bit worse for 
the non-smeared gauge background, where the gauge links at oblique angles systematically 
yield somewhat lower values. This error can be reduced considerably with our taxi driver 
correction, see section 4.4.8. Notice, however, that the amplitude decays much more quickly 
with increasing \£\ on the unsmeared ensemble than on the smeared one. The Gaussian-like 
shape emerging on the unsmeared lattice is much narrower than on the smeared one. We will 
address this issue in section 4.4, when we discuss renormalization. 
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(a) (b) 




Figure 4.2: Selection of link paths. 

(a) Illustration of the set of quark separations £ chosen for the analysis at P = 0. Only the 
gauge paths with £3 = are shown. 

(b) For the analysis at P = (—1,0,0) x 2ir/L, we choose the quark separations £ in such a 
way that the (£ 2 , ^-P)-plane is approximately covered evenly in the accessible region. 



4.3.2 The Nucleon at Non-Zero Lattice Momentum 

Let us also have a look at the results for FL fi (P, C^ at ) at P = (—1, 0, 0) x 2tt/L. Here we can 
calculate the amplitude A2 nien (£ 2 , £-P) at non-zero values of i-P. However, since we cannot 
implement link-paths with a Minkowski-time component £° 7^ 0, we are confined to the region 

£ 2 = -£ 2 < 0, \£-P\ < \P\ sf^P . (4.29) 

The above equation cuts out a wedge-shaped region in the (£ 2 , £-P)-plane, compare Fig- 
ure 4.2b. Only this region is accessible to us on the lattice. The wedge has an opening 
angle proportional to P, i.e., our investigations are limited by the highest nucleon momentum 
available to us. In order to be able to explore the accessible region thoroughly, we choose 
additional gauge links. Restricting ourselves to £2,-^3 > 0, we have determined a set of gauge 
paths such that the (£ 2 , £P) -plane is densely covered, see Fig. (4.2b). Performing the analysis 
for this extended set of over 1000 link paths on the HYP smeared ensemble results in the plots 
of Fig. 4.4. The dominant feature of the real part is the Gaussian-like decay with \/—£ 2 . The 
£-P-dependence is rather weak. The situation is different for the imaginary part. Here the am- 
plitude fulfills within statistics nicely the constraint Im A% nTen (£ 2 , £-P) = -Im A% nren (£ 2 , -l-P) 
following from eq. (4.6). The £-P-dependence is related to the rr-dependence of the TMD 
PDFs, and we will be concerned with it in detail in section 4.7. 
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(fm) 



Figure 4.3: Unrenormalized amplitude lA^^^l 2 , 0) obtained directly from the ratio 
Rj^,u-d(Pj C 1 ^), for P = on the coarse-m050 lattice. For clarity, we have averaged over link 
paths of the same shape (equivalent paths under H(4) transformations, see section 3.1.13). 
Link paths coinciding with the lattice axes are marked with a blue cross; the red error bars 
belong to link paths at oblique angles. 

(a) gauge links evaluated on the HYP smeared lattice, 

(b) gauge links evaluated without smearing. 
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Figure 4.4: The unrenormalized amplitude A^ en (l 2 , i-P) obtained directly from the ratio 
R~fe,u-d(T~i P, C\ at ) using the sequential propagators with P = (—1,0, 0) x 2ir/L on the coarse- 
m050 ensemble and applying HYP smearing to the gauge fields. Statistical error bars are 
shown as small squares floating above and below the interpolating surface, 
(a) real part, (b) imaginary part. 
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4.3.3 Restrictions from the Euclidean Lattice 

Figures 4.2b and 4.4 give us a vivid impression of the limited range of £ 2 and t-P accessible 
on the Euclidean lattice according to eq. 4.29. What are the respective implications for the 
calculation of TMD PDFs? From the formalism developed in section 4.1.3 we have learned 
that the Fourier transform of our amplitudes Ai(t,t-P) to momentum space always involves 
an integral over £ 2 from — oo to 0, and an integral with respect to t-P over the whole real 
axis, as in, e.g., eq. (4.17). Clearly, the lattice calculations cannot provide data in the whole 
integration range. 

Nevertheless, the amplitudes Ai(£,£-P) at t-P = give access to the first Mellin moments 
we have introduced in eq. (2.13). As we will show in more detail in section 4.6, we can 
carry out the Fourier transformation with respect to t 2 and get first Mellin moments such 
as /^(fcjj = ^[T + K 1 )(fej_; P 5 S). The first Mellin moment contains information about the 
probability densities of quarks with respect to transverse momentum k± alone. For example, 
in the case of /i we obtain the difference between the unpolarized quark and antiquark 

density. We can also learn something from the observable ^-independence of the amplitudes 
from the lattice, if we are willing to make additional assumptions, see section 4.7. 

However, without further assumptions, we cannot directly determine the x-dependence of 
quark densities. In particular, it is interesting to observe that PDFs, i.e., the fcj_-integrated 
distributions, are inaccessible to us. For example, from equation (4.17), we get 3 

d 2 k ± h(x,k ± ) " = h{x) = [ er i(i ' P)x 2A 2 (£ 2 ,£-P) o . (4.30) 



2tt 



e 2 =o 



Clearly, we have no freedom to vary t-P for £ 2 = on the lattice. Nevertheless, it is possible 
to calculate Mellin moments of PDFs, which corresponds to an expansion in terms of local 
operators. 



4.4 Renormalization of the Gauge Link 



As mentioned before, the differences between Fig. 4.3a and 4.3b indicate a strong renormal- 
ization scheme dependence of our operators. What are the renormalization properties of the 
Wilson line? 



4.4.1 The Wilson Line in the Continuum 

Let C be a continuously differentiable ("smooth"), non-overlapping contour of total length 
l[C]. In the continuum, the Wilson line along such a path is renormalized according to 



3 We remind the reader that, without special precautions, fi(x) is strictly speaking not the fcx-integral of 
fi(x,k±), see the discussion starting in section 2.3. 
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[DV80, Are80, CD81, Dor86] 

U rcn [C] = Z- 1 exp(-(5m/[C])^[C] . (4.31) 

Here Z z and 5m are renormalization constants. 4 The exponential factor in the equation 
above is an example of a power divergence, i.e., an ultraviolet divergence which behaves 
like a power of the renormalization scale. For example, for the regularization prescription 
1/x 2 — ► l/(x 2 + a 2 ), having an inherent renormalization scale ~ a -1 , Dorn [Dor86] finds at 
one loop order 5m oc g 2 ja + 0(g 4 ). In lattice gauge theory, which corresponds to a cutoff 
scheme with cutoff scale a -1 , perturbation theory gives a result of the same form, see section 
4.4.3 below. The power divergence also appears in other schemes, such as the Pauli-Villars 
scheme [PV49] and cutoff schemes in general. In dimensional regularization, 5m is zero, but 
renormalon ambiguities appear, compare, e.g., Ref. [Pin06]. 

Each "disruption" in the smooth Wilson line gives rise to another renormalization factor. 
Consider Wilson lines with a finite number n cusps [C] of cusps (i.e., points where the contour 
is not continuously differentiable) . Then the generalized formula for the renormalization of 
the Wilson line reads 



U Icn [C] = Z^exp -5ml[C] - V v(0i) U[C] , (4.32) 




where the renormalization constants v(Qj) depend on the opening angles 0{ of the cusps. On 
the other hand, for the color trace of a closed contour Ci oop , the multiplicative renormalization 
with Z~ l does not appear: 

(^cusps[Cloop] \ 
-5ml[C loop ]- Jl K^i)) tr c W[C loop ] . (4.33) 

For the operator Or°p(£) defined in eq. (4.1), renormalization factors arise at the end points 
of the Wilson line from quark field renormalization and the quark - gauge link vertex: 

0&(£) = Zfo g) Z- 1 ex V (-5ml[C}) O^l) . (4.34) 

V v ' 

For a straight gauge link, l[C] = \/—£ 2 , so the renormalized operator is of the form we have 
stated in eq. (4.27). It is of primary importance to get the length dependent renormalization 
with 5m under control, and we shall focus on this issue on the following pages. We will discuss 
the multiplicative renormalization factor Z^ z in section 4.6. 



4 This formula can be derived on very general grounds using the auxiliary z-field method [GN80] and 
BRS transformations [BRS75]. Wave function renormalization of the auxiliary field z is the origin of the 
renormalization constant Z^ 1 . 
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4.4.2 Renormalization Conditions 

The renormalization constant 5m cannot be determined unambiguously without an additional 
renormalization condition. To see this, let us write all a-dependences explicitly, and let us 
make the replacement 

5m(a) — ► 5m id(a) = 5m(a) + 5mo (4.35) 

in eq. (4.31). Then 

U Tea [C] = Z-\a) exp(-5m(a) l[C] - 5m l[C\) U[C]{a) . (4.36) 

The second /-dependent term is not a-dependent, so it can be absorbed in the definition of 
the renormalized Wilson line, and we get 

U ren2 [C] = Z-\a) exp(-5m(a)l[C}) U[C](a) (4.37) 

in a new renormalization scheme labelled "ren2" . Obviously 5m{a) is determined only up to 
a scale independent constant. First of all, we note that the dimensionless quantity Am(a) 
defined by 

Am(a) = a 2 —5m(a) = a— — 5m(a) (4.38) 
da a ma 

is free of the aforementioned ambiguity and thus adequate to specify results in a renormal- 
ization scheme independent way. For the comparison of ensembles with two different lattice 
spacings a\ and 02, it will be useful to approximate the derivative in the definition above by 
a finite difference. Here we choose to discretize the logarithmic derivative, i.e., the relative 
change in a: 

Am(ai, a 2 ) = ^~a~^ M^) - 8m(a{) _ ^ 

m(a 2 /ai) 

To fix a value for <5m(a), we need to specify a renormalization condition, i.e., we must provide 
some piece of information that uniquely defines the numerical values of the renormalized Wil- 
son line at lengths l[C] 3> a. In sections 4.4.7 and 4.4.5, we explore two such renormalization 
conditions. 



4.4.3 Perturbative Link Renormalization 

It is instructive to follow Refs. [EH90, BLP89, MMS92] and to use lattice perturbation theory 
at leading order to calculate the renormalization constant 5m for the discretized Wilson line 

Before we start, we need to introduce a lattice gluon propagator. To this end, the lattice 
gauge action S^ffJ] is expanded to lowest order in terms of A-fields using U(x,x + p) = 
exp(iga^(x)). Adding a gauge fixing term S^fC/] (necessary in perturbation theory) and 
expressing everything in momentum space, the action can be brought into the form 

S^[U) + S^[U) = l - J^^^Al{k)^(^b-^ (4.41) 
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(a) (b) 




X 



w T (j) 




Figure 4.5: Leading loop diagrams for gauge links in lattice perturbation theory, 
(a) Sunset diagram, (b) tadpole diagram. 



where kp = 2 sin(a£^/2) and where £ is the gauge fixing parameter. The expression for Dnp(k) 

,, . , _ ._ id nofinorl + 

'up — vab^vp 



is specific to the gluon action used. Now the gluon propagator Gp b ~ = 5 a bG v p is defined through 



the inversion 



Y. „2 (p~A k ) - - *A«*) G P p(k) = Sf^ (4.42) 
and the Fourier transformed propagator is defined by 

Gf u (x) = 5 ab j n/a ^expiikpxp) G^{k) = S ab £ ^^ e xp(^/a) G pp (p/a) , 

(4.43) 

where we have introduced the dimensionless momentum variable p = ak. Let us now expand 
the gauge link in terms of A-fields. For simplicity, we pick a path on one of the lattice axes, 
i.e., we set = nfi. We get 

n— 2 n— 1 -. n— 1 

(U^[C\) = 1 ~9 2 a 2 £ £ ^(*«H- t (*W) -i 5 V £ ^(,«) ^(z«) +0(/) 

i=0 j=i+l i=0 
s / v „ ' 



"sunset", see Fig. 4.5a "tadpole", see Fig. 4.5b 

^ n— 1 

= 1 - - 5 V ^ G%{x^ - *«) T a T fe + 0(g 4 ) , (4.44) 
i,i=o 

where T a , T b are the generators of SU(3) C . In terms of the gluon propagator in momentum 
space, we are able to perform the sums over i and j, and we obtain 

(u^n - 1 - i^ f ^-^h^y«)G-^m +0{gl) (4 45) 



(2vr) 4 1 - cos(p ( 

where l[C ut ] = an is the length of the gauge link and 1C F = ^ a T a T a = (4/3)1. 

Taking into account that G ^{p / a) / a 2 is o-in dependent, we take the derivative of the above 
expression with respect to a in order to isolate the divergent part: 



d / 7 ,latr~lat n _ tC F gH[C^\ r dp-, sm(p p l[C^]/a) [* d 3 p G m {p/a) 4 
daV [C ] >" 2a> ]_A2*Y l-cosfe) y_,(27r)3 ^ + °k ) ' 

(4.46) 
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S3 



where the three dimensional integration is carried out over all indices other than /i. The 
integrand of the p^-integral is a representation of the ^-functional, so that we get for a — ► 

Thus we find 

<«-(*-]> = <, - f £ (0^|^ o + n„n-d,ve rg e» t ) > + . 

(4.48) 

From this we identify the renormalization constant 8m: 

g 2 C F r d*p Gn-Jp/a 



8m 

a 2 



r d A p Gnn(p/a) 

/ To ^3 2 + const. 4.49 



= X 

Using eq. (4.38), this corresponds to Am = —g 2 X. Indeed, we find that 8m grows linearly 
with the cutoff a -1 . 

We can also take HYP smearing of the gauge fields into account. Following Ref. [DeG03], 
the effect of smearing can be implemented in our one-loop calculation by replacing the gauge 
field by a smeared one: 

A n {x) -> A^(x) = J2 KM M* + y) . (4.50) 

Here, the coefficients hpp(y) are specific to the smearing procedure. Effectively, this just 
means that we have to replace the gluon propagator eq. (4.49) according to 



G^(k) — > G^\k) = Y, hp'(~ k ) GM k ) hpip(-k) . (4.51) 

fi'p 



We have evaluated equation (4.49) numerically, both for the smeared and the unsmeared case, 
using the following ingredients, which we have listed for convenience in section B.2: 

• the inverse gluon propagator of the MILC action [Bis05] , 

• the parameters ci, C2 and C3 of the MILC action, where we have used the values for uq 
listed next to the unsmeared ensembles in Table (3.1), 

• the HYP smearing coefficients h^{k) [DeG03] with the values a% specified in section 
3.1.15, 

• for the coupling g, the bare lattice coupling listed in Table 3.1. 



The results are independent of the gauge fixing parameter £ and are listed in Table (3.1). As 
an interesting qualitative feature, notice that Am is much smaller on the smeared ensembles. 
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Note that the MILC gluon propagator reduces to the Wilson gluon propagator if we set a = 0. 
In this case, we get X = —0.168487, in agreement with Ref. [EH90]. 

Two improvements of the calculation discussed above are possible: Firstly, instead of using the 
bare lattice coupling, there are more sophisticated ways to adjust the coupling g. Secondly, we 
could go to higher loop order. Both of these strategies have been followed, e.g., in Ref. [MS99] . 
However, our aim in the sections to follow will be to examine non-perturbative methods to 
fix dm. What we have gained from the perturbative calculation is the confidence that our 
discretized operator will indeed need a length dependent renormalization of the same form as 
derived for the the continuum operator, and the insight that gluon self energy diagrams are 
responsible for the divergence proportional to the cutoff a -1 . 



4.4.4 Wilson Loops: A Study on Multiple Scales 

A Wilson loop 

W(r,t;a) = ^tv c U lat [C^} (4-52) 

is the color trace of a closed gauge path C!f| describing a rectangle of dimensions rxt (in phys- 
ical units) on the lattice. The major advantage of using Wilson loops for the determination 
of the renormalization constant 5m{a) is their gauge invariance. 

Can we determine the renormalization constants from the scaling behavior of the Wilson loop? 
According to the continuum formalism eq. (4.33) 

W ren {r, t) = exp (-tfm(a)Z - 4i/j_(a)) W(r, t; a) , (4.53) 

where we abbreviate 2(r + t) = I and the renormalization constant for 90° corners with 
v±_{a) = z/(90°; a). Demanding that W ren (r, t) be the same at two different lattice spacings a\ 
and a 2 , we get 

ln ( /w! r '! ,a2 !\\ J = *{M«g) - Mai)} +4 KM ~ "±(ai)} ■ (4.54) 

Obviously, we can determine the difference of renormalization constants at different lattice 
spacings. Let us define a quantity which can be expressed in terms of Am(ax, a 2 ) as introduced 
in eq. (4.39): 

v l * \ - 1 ( (W(r,t,a 2 )) \ 4yoT^ 

Y scsl {r,t;a 1 ,a 2 ) = jt— — ^rhi „ = Am{a 1 ,a 2 ) -\ : Au ± (a 1 ,a 2 ) , 

< ln(a 2 /oi) \(W(r,t,ax)) I I 



where Az/j_(ai, a 2 ) = {^j_(a 2 ) - v±_(a\)}/ ln(a 2 /ai) 



(4.55) 



To study the quantity Y scll i(r, t; oi, a 2 ) on the lattice, we have evaluated Wilson loops on 
ensembles which differ only in their lattice spacings. From the lattices listed in Table 3.1, 
we selected the supercoarse, coarse-m020, fine and superfine ensembles, which all have ap- 
proximately the same physical strange quark mass and a light to strange quark mass ratio 
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^u,d = 0.4m s . Using planar Wilson loops, only integer dimensions f x t are available on 
the lattice, making it difficult to compute loops of the same physical size for different lattice 
spacings. To overcome this, we interpolate In (W{r, t; a) J linearly in the (r, i)-plane. Typ- 
ical results for Y sca \(r, t; a\, 02) are shown in Fig. 4.6a and 4.6b. Clearly, for larger loops, 
^scai( r ; t] a i, °2) approaches a plateau value. 

We now evaluate Y sca x(r, t; ai, 02) on a grid of points in the (r, t)-plane. We restrict ourselves 
to the region r,t > 3max(ai,a2), where the linear interpolation of In (W^(r, t\ a)) works 
reasonably well. Also, we reject points with statistical errors of more than ±0.02. The 
results, plotted with respect to I = 2(r + £), are shown for the smeared ensembles in Fig. 4.6c. 
As they should, data points with the same I coincide, even if they differ in r and t. Using 
eq. (4.55), we can now determine Am(oi, 02) and Ai/j_(ai, 02) from fits to data from pairs of 
ensembles with similar lattice spacing. (In this case, we take uniform fit weights for all input 
data.) The results for Am(ai, 02) are given in the column labelled "scaling" in Table 4.1. We 
do not quote uncertainties. 

Notice that the numbers for Am(ai, 02) we get are much larger than the values Am(a) deter- 
mined perturbatively at similar lattice spacings. We conclude that our simple perturbative 
calculation gives very inaccurate results. However, some qualitative features agree: Am in- 
creases with the lattice spacing, and smearing reduces Am significantly. 

In the following sections, we will discuss methods that allow us to fix the renormalization 
constant 5m for a given ensemble. The numbers we have obtained for Am(ai,a2) from the 
scaling behavior of Wilson lines can serve as a valuable cross check of these methods. 



4.4.5 Renormalization with Wilson Lines 

Martinelli and Sachrajda [MS95, CGMS95] suggest to analyze Wilson lines along straight 
contours C\ of length I as follows 

, , 1 tr (lA^ [C lat 1 \ 

r -«> " -^(m) = -*»- ***.<r-m - - a v^ [c iT|) ■ <««> 

The expectation values of open Wilson lines are not gauge invariant; the equation above is 
therefore only valid in the context of gauge fixed ensembles. Here we choose Landau gauge 
fixing, see section 3.1.10. Our results for Yn ne (l) are plotted in Fig. 4.7a. We have selected 
data with errors of less than 0.02/a obtained with straight link paths on the axes of length 
I < L/2. According to the above equation, renormalization of the Wilson lines manifests itself 
as a shift of Yn ne (l) by a value 5m, which can be different for each ensemble. Thus we should 
shift the data points of each ensemble up or down until they all agree with each other. The 
agreement can be optimized at a certain value of I, i.e., at a certain renormalization point. 

For example, a suggestion in Ref. [MS95] is to adjust the curves at I — * 00, imposing the 
renormalization condition lim^oo 4 In tr c U Ten [C{\ = 0. Note that the authors point out that 
the existence of this limit is theoretically not proven. Fitting the form aln ne (/) = —5m — 
7m(l + a/(l — a/2)) (as suggested in Ref. [CGMS95]) to data with I > 5a, we arrive at 
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(a) (b) 




o.oo 

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 
2(r+t) (fm) 



Figure 4.6: (a) and (b): Y sca \(r, t; oi, 02) evaluated from the ratio of equally sized Wilson loops 
on the coarse-m020 and fine ensemble. We only include points with errors below ±0.02. (a): 
with HYP smearing, (b) without HYP smearing; statistics breaks down much earlier, 
(c): Y sca i(r, t; ai, 02) from pairs of HYP-smeared ensembles at m u ^ = 0.4m s . Selected data 
points fulfill r, t > 3a (a taken from the larger lattice) and have statistical errors below ±0.02. 
The fits are performed with the parametrization given in eq. (4.55). 
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Fig. 4.7b. In this plot, we have already offset the data by 5m, so that the plotted fit functions 
approach zero at infinity. The fit is not very stable, and the values Am we obtain from this 
method are not in good agreement with our scaling analysis in the previous section. 

Therefore, let us try another renormalization point, and demand 4 lntr c £/ ren [C/] = at 
I = 0.5 fm, a scale at which we have accurate data. Using simple linear interpolation between 
data points at ^-values above and below 0.5 fm, we can determine 5m and arrive at Fig. 4.7c. 
Numbers are listed in the column labelled "Wilson line at I = 0.5 fm" in Table 4.1. We have 
also listed the resulting values Am{a\,a2) between adjacent ensembles, and find very good 
agreement with the scaling analysis of the previous section. 

In Fig. 4.7d, we show Yn ne (/) with offsets 5m determined from renormalization with the string 
potential, as described in the following section. In our case, this method obviously does not 
perform very well on the unsmeared ensembles: As we shall see, this problem is not unexpected 
and could be resolved with more input configurations. Therefore, we restrict ourselves to the 
smeared ensembles in Fig 4.7e. 

Looking, in particular, at Fig. 4.7c and Fig. 4.7e, we find that the data of different ensembles 
can be brought to almost perfect agreement for Wilson lines of lengths greater than about 
0.3 fm. This observation corroborates our notion that the corresponding lattice operators 
approximate continuum operators with straight Wilson lines, and thus share their renormal- 
ization properties. On the other hand, it is obvious that our data exhibits substantial lattice 
cuttoff effects for Wilson lines of lengths < 0.25 fm. For such short Wilson lines the continuum 
inspired renormalization prescription fails: the data points of different ensembles do not just 
differ by a constant shift 5m. Hence we will discard data for operators with links shorter than 
0.25 fm in our analysis of nucleon observables. 

Moreover, we learn that is important to impose a renormalization prescription which is sen- 
sitive to the data in the region where statistical and systematic errors are both small: The 
gauge link should be long enough, so that finite-a-effects are small, and yet not too long, so 
that finite volume effects and statistical uncertainties are acceptable. 

Finally, we draw again attention to the fact that different renormalization conditions may 
produce values for 5m that can differ by a renormalization scale independent constant. This 
is the reason why the plots in Fig. 4.7b-4.7d feature different offsets in the ordinates. 

4.4.6 Discretization Errors estimated with Wilson Lines 

A hint about the size of lattice cutoff effects can be obtained from the comparison of smeared 
and unsmeared ensembles. To this end, let us take a look at 

A[5m] = [Yc e (h) - rcrm - [Yc e (h) - rcrrn . (4.57) 

In each of the square brackets we calculate the difference between the renormalization con- 
stants 5m on the smeared and unsmeared ensembles, albeit at two different renormalization 
points l\ and fo- If we had no lattice cutoff effects, we should obtain exactly the same result at 
both renormalization points, and A [5m] would be zero. To be able to evaluate Y^^l) for any 
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Figure 4.7: The length derivative of the straight Wilson line in Landau gauge Y]i ne (l) from the 
four ensembles at rh U( i = 0.4m s , with HYP smearing (filled symbols) and without smearing 
(open symbols), (a) unaltered data, (b) fitted and renormalized based on an assumption 
about the asymptotic behavior, (c) renormalized by imposing a renormalization condition at 
I = 0.5 fm, (d) renormalization based on the string potential, see section 4.4.7, (e) same as the 
previous item, but for the smeared ensembles only. The data in the shaded region is clearly 
affected by strong lattice cutoff effects. 
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I on a given ensemble, we generate spline interpolations. Setting l\ = 0.25 fm and I2 = 0.75 fm, 
we obtain A[5m] = 0.027(11) GeV, A[5m] = 0.0308(26) GeV, A[Sm] = 0.0384(14) GeV, 
A [5m] = 0.04711(56) GeV on the superfine, fine, coarse-m020 and supercoarse ensembles, 
respectively. These numbers suggest that A[5m] extrapolates linearly to at a = 0. 

For our calculations of nucleon structure in later sections, we will work with data obtained 
on the coarse ensemble with quark separations ranging from about 0.25 fm to 1.2 fm. So 
let us set h = 0.25 fm and I2 = 1.2 fm. Within our limited statistics, we cannot determine 
A [5m] with reasonable precision on the two finest lattices for these choices of l\ and li- 
However, on the coarse-m020 ensemble, we obtain A [5m] = 0.0393(88) GeV, or, in lattice 
units, A [5m] = 0.0238(54). This number will serve as an estimate of systematic uncertainties 
from lattice cutoff effects. We also calculate this error on the coarse-m050 ensemble, where 
we obtain A [5m] = 0.0405(56). The uncertainties thus obtained should not be ignored; they 
are in general much larger than uncertainties from other sources. 

4.4.7 Renormalization based on the Static Quark Potential 
4.4.7.1 The Principle 

The two Wilson lines running parallely in t direction in a Wilson loop can be interpreted 
as propagators of static quarks Q and Q, compare section 5.1 and textbooks such as Refs. 
[Rot97, DeG03]. The other two Wilson lines are then regarded as a gauge invariant source of 
a QQ pair separated by a distance r. Just as we did it for the nucleon in section 3.2.4, we 
can extract the energy of the QQ system from the slope of the logarithmic correlator at large 
t. We call this energy static quark potential 

V(r; a) = - lim — In (W(r, t; a) » . (4.58) 

Replacing V(r;a) and W(r,t;a) by renormalized quantities, and substituting eq. (4.53), we 
obtain 

d 

V ren (r) = 2 Sm(a) - lim — In (W(r, t; a)) = 2 5m(a) + V(r; a) , (4.59) 
t^oo at 

= 2m r Q en -£ bind (r) , (4.60) 

which has an interpretation in terms of a binding energy -EWd( r ) and the mass of the static 
quark rriq 11 in the last line of the above expression. We observe that a change in the renor- 
malization condition for the Wilson line 5m(a) — > 5m(a) + 5mo is equivalent to a finite quark 
mass renormalization m™ — > m™ + 5mo/2. 

We can fix 5m(a) by imposing a renormalization condition on V ren {r). We want the condition 
to be insensitive to our determination of the lattice scale a. One successful approach [C + 08, 
Pet07, Baz08, B + 09] makes use of the observation that the static potential at large distances 
r can be very well described by the formula 

^frfngM = or - ^ + (4.61) 
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derived in Ref. [LSW80]. The formula 5 originates from a small-h approximation based on 
the assumption that the Wilson loop satisfies the equation of motion of a quantized string 
[Nam79]. Note that an exact result for the string potential is also available [Arv83], which 
shows that the string model certainly cannot describe the static quark potential for smaller 
values of r. In the expression above, the string tension a is a fundamental constant, while 
C ren is the renormalization constant related to the self energy of the Wilson loop. It turns 
out that the lattice data can be fitted very well (in terms of dimensionless quantities) to a 
function of the form 

aV(r; a) « V &t (r; a) = & a r - °^ + C a (4.62) 

r 

with the fit parameters a a , a a and C a . Let us set the lattice potential and the string potential 
equal at some fixed r = r matc h, which we choose large but still in the range where lattice data 
is available. Taking a = cr a /a 2 and putting together eqns. (4.60), (4.61) and (4.62), we get 

VfiHrWdO = V; r tring(^match) => 2 8m{a) = -°— (a a - £-) + (aC ren - C a ) . (4.63) 

^match ^ -I ^ ' 

Here the lattice spacing a appears explicitly in combination with r matc h, which we must specify 
in physical units in order to implement the renormalization condition scale-independently. It 
is useful to eliminate this appearance of a in favor of the Sommer scale ro, which is defined 
by the condition [Som94] 

2 dV(r;a 



Oi 



= 1.65 . (4.64) 



Effective potentials reproducing the experimental spectrum of heavy quarkonia (bb bound 
states, etc.) show that ro ~ 0.5 fm. 6 Using our parametrization eq. (4.62), the Sommer scale 
in lattice units is f-Q = a/ (1.65 — a a )/a a . Substituting a = tq/tq we finally obtain 



2 Sm(a) = ( ) J = ° a («. - - C a + aC™ . (4.65) 

V r match/ V 1-65 - a a V 12/ 

In order to be able to use this equation to determine 2 8m(a), we need to make a choice 
regarding C ren . This choice is part of the definition of our renormalization condition and 
is related to the arbitrary shift 5mo in eq. (4.36). We will follow Ref. [C + 08] and choose 
C ren = 0. The right hand side is now free of any explicit a-dependence. Together with the 
convention r matc ^/rQ = 1.5 of Ref. [C + 08], the equation above forms a practical and robust 
renormalization prescription. 

The approach discussed above may be interpreted in the following way: Matching to the string 
potential, we anticipate that the static potential will converge to a straight line or + C ren 
for large r. 7 By setting C ren = 0, we demand that this line run through the origin, see the 
dashed straight line or in Fig. 4.9a. 



The excellent agreement of KtrfngM with lattice data at large r is shown in Fig. 1 of Ref. [C 08]. 
6 Note that the MILC collaboration prefers to modify the condition above, replacing the constant 1.65 by 
1. The corresponding scale ri = 0.317(10) fm [A + 04] has been used to determine the lattice spacings quoted 
in Table 3.1. 

7 Just from the symmetry of the Wilson loop under exchange of temporal and spatial extent, it can be shown 
that the Wilson loop is bounded from above by a linear function in r [Sei78]. 
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Figure 4.8: Rectangular Wilson loop in the calculation of the static quark potential. 

(a) Wilson loop with step-like, spatial connections at an oblique angle. 

(b) Gluon exchange between the temporal Wilson lines in leading order lattice perturbation 
theory. 



4.4.7.2 Implementation and Results 



Apart from planar Wilson loops, let us consider Wilson loops with oblique, step like gauge 
paths in the spatial direction, see Fig. 4.8a. The temporal Wilson lines of these loops will then 
be separated by a lattice vector r in the spatial hyperplane. To calculate the static quark 
potential V(r), we pick a vector r and select Wilson loops with a temporal extent that is at 
least t m i n (see Table 4.1) and at most T/2. Assuming that the ground state dominates for 
these temporal extents, we can perform a fit of the form 



(W(r, t; a)) = b r exp (-iv r ) , 



(4.66) 



where b r and v r are fit parameters. We set the fit weights to the bootstrap errors of the 
individual input data points. 

The values v r obtained in the fits correspond to the potential at the given r = \r\. At small |r|, 
the potential thus determined exhibits considerable breaking of rotational invariance, which is 
a clear sign of discretization errors. However, these errors can be considerably reduced with a 
perturbative correction [Mic92, BalOl]. Calculating the gluon exchange between the temporal 
Wilson line segments in the loop (cf. Fig 4.8b) to leading order in lattice perturbation theory 
similarly as in section 4.4.3, one obtains 



lim — 

t^oo dt 



In (W(r, t; a) 



C F 9 2 



d 3 p . G u (p/a) 
3 exp^r^) - 2 



CfI 
4tt 



(4.67) 



The corresponding calculation in continuum perturbation theory simply produces the Coulomb 
potential Cpg 2 / '(4-7r|f |). Indeed, [1/r] approaches 1/f with increasing r very quickly, see ap- 
pendix B.3 for details. To apply the correction in practice, choosing the appropriate coupling 
g is non-trivial, so the usual method is to determine a strength parameter A a for the lattice 
artefacts along with the parameters a a , C a and a a in a fit to the potential. In this fit, the 
values v r determined above serve as input data. The fit constraits are of the form 



OL 



cr n r 



+ C a — X a 



1 



1 



(4.68) 
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Again, we choose fit weights according to the statistical errors of the input data. Taking the 
fit results thus obtained, we can determine the renormalization constant 5m{a) according to 
eq. (4.65). The parameter A a does not enter the renormalization condition. However, it turns 
out that the additional degree of freedom introduced through the perturbative correction is 
very important for the quality of the fit eq. (4.68). Only with the corrections, we can use 
data points v r down to r = y/2. For details regarding the corrections, see section B.3 in the 
appendix. As a consistency check, we also calculate the lattice scale a determined from the fits. 
On the smeared lattices, we obtain deviations below 3% with respect to the numbers obtained 
from the MILC collaboration. Systematic errors could be further reduced by increasing the 
minimal size of the Wilson loops, which would however demand a larger number of input 
configurations. We have estimated these systematic errors for 5m by varying the minimal size 
of the Wilson loops in the fits and get errors that are negligible compared to the systematic 
uncertainties determined from the smeared-unsmeared comparison in section 4.4.6. The final 
results for 5m are quoted in Table 4.1. 

Figure 4.9a shows the static quark potential after renormalization with the method described 
above on the HYP smeared ensembles. The data points of the different ensembles lie very 
close to each other. This is a clear indication that the method works: The renormalized static 
quark potential from the lattice is unique; in particular it is independent of the lattice spacing. 
Without renormalization, we would see large offsets between the lattice results for the static 
quark potential on the different ensembles. 

Our results are worse on the unsmeared ensembles, see Fig. 4.9b. However, it is obvious 
that our calculations for the unsmeared ensembles could be easily improved by taking more 
configurations (to improve statistics) and by raising the minimal temporal extent of the Wilson 
loops t m i n (to reduce systematic uncertainties). 



4.4.8 Taxi Driver Correction 

In the previous sections, we have discussed ways to determine a constant 5m which can be used 
to renormalize straight Wilson lines on the lattice. For a straight link path Ci at connecting 
two points separated by a vector £, the length is related to the number of link variables by 
\£\ = l[C l&t ] = a ni; n ks[C lat ]. What about the step-like paths that we use to discretize Wilson 
lines at oblique angles? We have already observed in Fig. 4.3 that operators with step-like 
link paths produce expectation values that lie systematically lower than those of operators 
with straight link paths of the same length \£\, especially on the unsmeared ensembles. Like a 
taxi driver navigating on a rectangular grid of streets and avenues, the link path on the lattice 
has a longer length than the continuum Wilson line we wanted to model: ani; n k s [C lat ] > \£\. 
The link path has the shape of a polygon, which approximates a straight, smooth Wilson 
line and has an infinite density of kinks in the limit a — > 0. In the continuum, Wilson lines 
of this shape have already been discussed in Ref. [CD81]. They are renormalized like the 
smooth Wilson line they approximate, but the renormalization constant 5m receives an extra 
contribution. Inspired by this continuum result, we may conjecture that the step-like Wilson 



94 



TMD PDFs with Straight Links 



(a) 




0.0 0.2 0.4 0.6 0.8 1.0 

r(fin) 



(b) 

1.0 

0.5 

> 

S o.o 

^ -0.5 
-1.0 

0.0 0.2 0.4 0.6 0.8 1.0 

r(fm) 

Figure 4.9: Renormalization using the string potential. The data points shown have been 
shifted by the renormalization constants determined from the fits and have been corrected for 
known discretization errors by subtracting A a ( [1 jr ] — 1 jv ) . The thin solid lines through the 
data points correspond to fits of the form eq. (4.62). For the shown string potential V^rfng 
from eq. (4.61), we take C ren = and an average a determined from the fits to the smeared 
ensembles. The vertical dashed line marks approximately 1.5ro, where the matching to the 
string has been performed. 

(a) HYP smeared ensembles. We also show the linear part ar of the potential for an averaged 
numerical value of a. 

(b) Unsmeared ensembles. The statistics we have accumulated here is very limited; the small 
values of i m i n we have used to get some coarse results may entail significant systematic un- 
certainties. For the red dashed line showing the string potential, we took the same parameter 
values as in the smeared case. 
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line on the lattice can be improved by a correction factor: 



W^ rr [C] at ] = exp (-(n links [C lat ] - \£\/a) 5m taxi - n cusps [C lat ] iw) Z/ at [C] at ] . (4.69) 



Here n cusps [C ] is the number of 90° angles in the link path C , and 5mt ax i and ft ax i are 
constants, which we simply determine from a fit, based on the requirement that the corrected 
operators should give expectation values that depend smoothly on \£\. To this end, we take 
expectation values of Wilson lines 



calculated on a Landau gauge fixed ensemble. We interpolate the results for straight link 
paths by a natural spline, as shown in Fig. 4.10a and 4.10c. Next, we apply the correction 
eq. (4.69) to the Wilson lines with step-like link paths, adjusting 5mt ax i and ft ax i until the 
mean squared distance between these data points and the spline becomes minimal. The fit 
weights are chosen according to the bootstrap errors of the individual data points. In the 
fit, we restrict ourselves to data points with 2a < \£\ < L/2 in order to exclude the regions 
where strong lattice cutoff effects and finite volume effects may be expected. The fit results 
are listed in Table 4.2. 

The corrected data points, shown in Fig. 4.10b and 4.10d, do not lie perfectly on the spline 
within their tiny errors. Obviously, the prescription eq. (4.69) cannot remove all artefacts 
created by the use of step-like link paths. Nevertheless, the improvement is significant, in 
particular on the unsmeared ensembles. For example, the weighted root mean square distance 
to the spline (the square root of "x 2 ") is reduced by a factor of almost 14 on the unsmeared 
coarse-m020 ensemble. The correction seems to work well for the whole range of \£\. Even 
the two data points in the region excluded from the fit in Fig. 4.10a are visibly improved in 
Fig. 4.10b. On the smeared ensembles, where the data points lie close to the spline from the 
start, the improvement is much smaller. 

We conjectured in Ref. [MHS + 08] that the constant <5mt ax i determined from the fit to the 
spline might be used for the overall renormalization of the Wilson line, i.e., that the entire 
renormalization of the lattice operator depends on the number of link variables and corners 
only: 



The numerical results of our studies at several lattice spacings are not in support of this 
conjecture. In particular, the values Arh(ai,a2) obtained from <5rn. tax i disagree with our 
results from section 4.4.4. They are too small by a factor of more than 2 on the unsmeared 
ensembles, and by roughly an order of magnitude on the smeared ensembles. It appears that 
the taxi driver correction merely addresses a certain kind of discretization error. On top of the 
taxi driver correction, an |£|-dependent renormalization as discussed in the previous sections 
is mandatory. 




(4.70) 




(4.71) 
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Figure 4.10: Taxi driver correction: We plot the unrenormalized expectation value of the 
Wilson line ^-Retr c (W lat [C lat ]} on the coarse-m020 ensemble. The blue crosses indicate 
results from straight link paths, the black dots are obtained with step-like link paths. Error 
bars are smaller than the symbols in this plot, 
(a), (b): without smearing, 

(c), (d): evaluated on the HYP smeared ensemble. 

In (a) and (c), we show the uncorrected data, and the smooth curve is a natural spline 
interpolating the blue crosses. The shaded area is the region excluded from the fit that 
determines (5m t axi and 5f tax j. In (b) and (d) we plot the data after the taxi driver correction. 
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ensemble 




^taxi 


super coarse smeared 
coarse-m020 smeared 
fine smeared 
super fine smeared 


-0.0300(12) 

n noo9 / 1 o\ 
-U.0zzo(izJ 

-0.02217(93) 

-0.02180(91) 


0.00728(45) 
0.00624(52 J 
0.00735(57) 
0.00762(53) 


super coarse 
coarse-m020 
fine 

super fine 


-0.2294(15) 
-0.20664(72) 
-0.18511(61) 
-0.1804(13) 


0.04700(65) 
0.04756(36) 
0.04387(31) 
0.04682(76) 



Table 4.2: Parameters determined for the taxi driver correction. 



4.5 Dividing Amplitudes by the Wilson Line 

Consider the operator 

n div m = Prop (i) = g(i)r op WA g(Q) u 79 , 

FOP[ j " ^Retr c (U[0J]} &Retr c <W[0,<|) ' 

where the expectation value of the Wilson line in the denominator requires some gauge fixing. 
We choose the Landau gauge. According to eqns. (4.31) and (4.34), the operator defined 
above is renormalized according to 

OgTM = z-: 1 z 2 m 0«%(£) . (4.73) 

The renormalization factor exp(— Sml[C]) cancels. Thus the unrenormalized amplitudes ex- 
plored in section 4.3 divided by the Wilson line expectation value 

Af v (£ 2 ,£-P) = -p V ' (4.74) 

are quantities that need no ^-dependent renormalization. In section 2.4.5 the square root of 
a Wilson loop serves as a subtraction factor with the same function as the Wilson line here. 
The obvious disadvantage of using the Wilson line as subtraction factor is that we introduce 
dependence on the gauge fixing condition. 

Figure 4.11 shows some examples of amplitudes divided by Wilson lines, plotted as a func- 
tion of \/—£ 2 . The amplitudes have been obtained from ratios as described in section 4.1.2. 
Surprisingly, the data shows only weak dependence on \J—£?. It appears that the charac- 
teristic Gaussian-like decay we saw in the unrenormalized amplitudes, e.g., in Fig. 4.3, has 
been cancelled in the ratio with the Wilson line. Some of the divided amplitudes are almost 
constant within the whole range of £ 2 accessible to us. These results suggest that the typical 
Gaussian-like shape of the unrenormalized amplitudes is not a feature characteristic of the 
nucleon. Instead, the Gaussian decay seems to be largely driven by the Wilson line, which 
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already shows this behavior when evaluated with Landau gauge fixing in the vacuum. This 
indicates that the ^-dependence of our straight link correlator is predominantly given by the 
dynamics of the gauge field. 



4.6 The First Mellin Moment 

In section 4.3.3 we showed that we cannot calculate the amplitudes for all (£-P)-vahies that 
enter the Fourier transformation to momentum space. However, if we restrict ourselves to 
the first Mellin moment as defined in section 2.3.3, the Fourier transform with respect to I 
is restricted to £ + = £~ = 0, i.e. we only need the amplitudes for t 2 < and t • P = 0. 
We plot the real parts of the amplitudes A2, Aq and Aj at l-P = evaluated on the HYP 
smeared coarse-m020 ensemble in Fig. 4.12. (The imaginary part of these amplitudes vanishes 
at l-P = due to eq. (4.6).) In the following, we explain how we arrive at the renormalized 
amplitudes and the Gaussian parametrization shown in this figure. 



4.6.1 Gaussian Parametrization 

In Ref. [MHS + 07], we showed that a type of Gaussian fit function can successfully describe 
our results for the unrenormalized amplitude A2(l 2 ,0) from the lattice (compare Fig. 4.3a) 
and enabled us to make the Fourier transform to momentum space. The Gaussian ansatz 
has been widely used and has a number of virtues [C + 06b]. A look at first Mellin moments 
obtained from models also suggest that a Gaussian function might be a good starting point, 
compare Fig. 2.10. Now that we know more about the renormalization of our operators, we 
still find the Gaussian parametrization useful, although there are certain limitations, which 
we shall address in the section 4.6.3. The simple Gaussian fit function we presently use for 
our amplitudes reads 

1,(^,0) « ^yexp = ^exp (-^j , (4.75) 

where Cj and aj are the fit parameters, the latter one characterizing the width of the amplitude: 



1/2 ^ [MUiA jt ei^ = 

V J Ml A 2 (-e 2 ± ,o) 3 v ; 



'A, 



Let us now study TMD PDFs based on the Gaussian parametrization. In terms of the 
amplitudes Aj, the first Mellin moment defined in eq. (2.13) is of the form 

*™ (fe±; P) S) = ±[0- 2 e*^ £ a^l^ e ±ni • • • i ±nmj 2A j{ f,0) , (4.77) 

j,ni,...,n mj 

where the a - m „ m are appropriate coefficients that can be read off from the parametrization 
eq. (4.5). Inserting the Gaussian model function eq. (4.75), we can readily perform the Fourier 
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Figure 4.11: The divided Amplitudes Af lv (£ 2 ,£-P = 0) evaluated on the HYP smeared en- 
sembles. In some of the plots, red horizontal lines have been inserted at an average value 
to guide the eye. Note that statistical errors are smaller on the coarse-m050 ensemble. The 
amplitude Af has been calculated only for the coarse-m020 ensemble up to now. (The shaded 
region indicates that we expect significant lattice cutoff effects in the amplitudes below about 
0.25 fm.) 
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transform and get 

(4.78) 

By comparison with eqns. (2.49) - (2.51), we can easily identify the first Mellin moments of 
TMD PDFs. In particular, we find 



47T \ 



( 2 M) 2 , 

= /i (1) (fc±), (4.79) 

*™« (Jfc±; P, 5) = -A^I exp ( -^-) - *^ exp ( 

4vr Vl 2 /^) 2 / raw 8vr Vl 2 /^) 2 



A 9 mki) + ^-^ 9 ^(ki) . (4.80) 
m N 



4.6.2 Renormalized Amplitudes from Gaussian Fits 

According to eq. (4.34), the renormalized amplitudes are given by 

Aj (£ 2 , 0) = Z^ z exp [-5m y/^Fj I™(^ 2 , 0) . (4.81) 

In a first step, we renormalize our data with respect to the £ 2 -dependent term, using the 
value 8m obtained with the help of the static potential, as described in section 4.4.7. We 
also include the taxi driver correction from section 4.4.8, although the effect is small on 
the smeared ensembles. Next, we fit Gaussians of the form eq. (4.75) to the data, from 
which we obtain the fit parameters Cj nren and Oj. The coefficients c" nren do not include 
multiplicative renormalization with Z7 z , yet. For the Gaussian fit, we exclude data points 
with \f— £ 2 < 0.25 fm, because gauge links of such short lengths are subject to significant 
lattice cutoff effects, which is an observation we made looking at Fig. 4.7e. We also exclude 
data points with \/—£ 2 > L/2, i.e., for operators larger than half the box size of the periodic 
lattice. For such operators, the separation between "copies" on the periodic lattice can become 
smaller than the extent of the non-local operators themselves. This might introduce finite 
volume effects. However, we remark that the results would not change much if we were to 
include the data points with \J — £ 2 > L/2 in the fits. In the fit procedure, we adjust the fit 
weights to the statistical errors of the individual data points. 

In the last step, we address the multiplicative renormalization factor z . Here we explicitly 
make use of the assumption that fi q (x, k±) has an interpretation as the density of quarks of 
flavor q in the proton, as discussed in 2.5.1. From this assumption follows 

J d 2 k ± j dx h, q (x,k ± ) = { I ; q q Z U d , (4.82) 
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where the right hand side specifies the number of valence quarks in the proton. Sea quark 
contributions cancel in the integration over all x, because the antiquark densities are given 
by fi,q{x, kj_) = —fi t q(—x, k±) [TM95]. For the determination of Zl z , we take the results for 
q = u — d, where disconnected diagrams in our three-point function cancel. Thus, in terms of 
amplitudes, we simply demand 

2i 2 ,«-d(0,0) = l . (4.83) 
In practice, we determine ZT 1 ^ and the renormalized coefficients cj according to 

rr „unren . ~— 1 „unren (a qa\ 

■- c 2,u-d> c j ■- Z ^, z C j ■ \ A - M ) 

Note that we deliberately do not use the data point at £ 2 = for renormalization. As 
mentioned before, we exclude it from the fit, along with all other results obtained with gauge 
links shorter than about two lattice spacings. The corresponding operators must be regarded 
as local operators; their renormalization properties are different and can be treated with other 
techniques than the ones discussed in this work, see, e.g., Ref. [MPS + 95, G + 99]. For local 
operators, the multiplicative renormalization constants depend on r op . In contrast, in the 
non-local case we consider here, the constant Zu z is universally applicable regardless of r op , 
see section 4.2 of Ref. [Dor86]. 

Numerical results for the Gaussian fits are given in Tables 4.3, 4.4a and 4.4b. In Table 4.3, 
the values for C2 show that the condition eq. (4.82) is approximately fulfilled in the u and d 
channel using the renormalization constant Zl z determined from the u — d channel. This is 
an indication that contributions of disconnected diagrams are negligible at the present level 
of accuracy for the observables we study here. 

4.6.3 A Critical Look at the Renormalization and Fit Prescription 

A serious issue in the above procedure is the choice of the renormalization condition for the 
Wilson line. In Fig. 4.12 we have renormalized with the static potential and C ren = as de- 
scribed in section 4.4.7. The picture changes drastically if we choose another renormalization 
condition. Just to show the effect, we plot in Fig. 4.13 data which has been renormalized 
with 8m obtained from the Wilson lines in a Landau gauge fixed ensemble, as described in 
4.4.5. This corresponds to the choice C mn ~ 0.2 GeV in terms of renormalization with the 
static potential. The renormalized data now keeps rising out to ^/—i 2 ~ 0.6 fm. Clearly, a 
Gaussian fit would not work in this case, unless we are willing to exclude a much larger range 
of data points from the fit. Is there a physical interpretation of the renormalization condition, 
i.e. the value C ren ? Can we establish a relation to a continuum factorization scale, or to a 
continuum renormalization scheme and scale? These are important questions that have to be 
addressed in future studies. 

Another issue concerns the high-fcj_-behavior of TMD PDFs we already encountered in sec- 
tions 2.3.2 and 2.5.3. The Fourier transformed amplitude, f[ (k±), is in turn a Gaussian 
function, see eq. (4.79). Obviously, a Gaussion does not adequately describe the large-fej_- 
behavior ~ predicted by perturbation theory. On the other hand, the fcj_-integral in 

our normalization equation eq. (4.82) is now well-defined, even without a \k±\ -cutoff. With the 
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Figure 4.13: Amplitude A2(£ 2 ,£-P = 0) evaluated with straight Wilson lines on the HYP 
smeared ensemble coarse-m020, with a different renormalization condition as in Fig. 4.12: We 
take 5rh from our renormalization procedure with Wilson lines on the Landau gauge fixed 
ensemble at a renormalization point I = 0.5 fm as described in section 4.4.5. The data is not 
renormalized with respect to the multiplicative factor Z^ l z . 



Gaussian approach, no explicit cutoff is needed. Instead, the Gaussian ansatz itself functions 
as a regularization prescription. 

At the level of amplitudes, the asymptotic behavior f[ (k±) ~ b/k\ translates into a con- 
tribution of the form —2nb ln(\/—£ 2 /lo) to A2(£ 2 ,0) at small \/—£ 2 , see section B.4. This 
contribution diverges at I = 0. Note that our lattice results do not exhibit the divergence 
at i = because the lattice imposes a momentum cutoff. Instead, we observe effects of the 
lattice cutoff, which is the reason why we exclude data points with \f-P- < 0.25 fm. Our 
Gaussian fit constitutes an interpolation that smoothly bridges the excluded region at small 
\[— ~£? . Roughly speaking, with our fit prescription, we "do not resolve" short range behavior 
at \J — £ 2 ;$ 0.25 fm. Correspondingly, we expect that the TMD PDFs we calculate become 
unreliable for k ± > 1/0.25 fm » 0.8 GeV. 

In the future, we should investigate whether we can detect, isolate and interpret the onset of 
the logarithmic short-range behavior in our amplitudes from the lattice. Special attention will 
have to be payed to lattice cutoff effects, which also set in at small \J — £ 2 . Fit functions that 
exhibit the correct behavior in the perturbative regime will enable us to treat the high-fcj_ 
behavior in a more systematic way. 



4.6.4 Comparing a Smeared and an Unsmeared Ensemble 

In Fig. 4.3 we observed a big difference in the widths of the unrenormalized amplitudes ^ nrcn 
on the smeared and unsmeared coarse-m050 ensemble. Does the renormalization procedure 
outlined here eliminate these differences? Up to now, we do not have renormalization con- 
stants from the static quark potential available for the coarse-m050 ensemble, which is the 
only ensemble for which we have calculated correlators with unsmeared gauge links. To be 
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Figure 4.14: Amplitude A<i{fi,l-P = 0) evaluated with straight Wilson lines for u—d quarks 
on the ensemble coarse-m050. We compare the results for gauge links evaluated on the 
unsmeared gauge configurations to the results obtained with HYP smearing. We take 8m 
from our renormalization procedure with Wilson lines on the Landau gauge fixed ensemble 
at a renormalization point I = 0.5 fm as described in section 4.4.5 and add 0.1945 GeV to 
reproduce approximately the renormalization condition from the static quark potential with 
C ren = 0. 

able to make the comparison nonetheless, let us use 5m as determined from the Wilson line 
at the renormalization point 0.5 fm and add a constant 0.1945 GeV. This will approximately 
correspond to renormalization with the static potential at C ren = 0. We show the results 
in Fig. 4.14. Smeared and unsmeared results agree much better, but there is still an ob- 
servable difference. The inverse widths of the Gaussian fits to the renormalized data are 
4/0-2 _ o.4146(85) sta t GeV for the unsmeared, and 4/cr| = 0.3872(34) stat GeV for the smeared 
ensemble. However, the discrepancy is no larger than we would expect it to be from our study 
of Wilson lines in section 4.4.6: With the help of eq. (B.18), we obtain a systematic error in 
4/o| of ±0.037 GeV. 

4.6.5 TMD PDFs and Densities from the Gaussian Parametrization 

Notwithstanding the issues of section 4.6.3, we give an interpretation of our Gaussian fit 
results in terms of TMD PDFs and quark densities. From eqns. (4.79) and (4.80) we readily 
read off and g^, in terms of the fit parameters and plot them in Fig. 4.15. In our 

plots, we label our profile functions with an extra superscript "sW", to remind the reader 
of the fact that these TMD PDFs have been obtained with straight Wilson lines and are 
therefore not strictly identical to the TMD PDFs defined and used in the literature and for 
the description of, e.g., SIDIS. 

As discussed in 2.3.3, our first Mellin moments are combinations of quark and antiquark 
densities. Presumably, the contribution from antiquarks is small, so it is interesting to test 
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Figure 4.16: A numerical study of bounds on TMD PDFs: On the HYP-smeared coarse-m020 
ensemble, we test whether the function A 9 (|fcj_|) defined in the text is positive. The error 
bands include statistical errors only. Renormalization is based on the static quark potential 
with the renormalization condition C ren = 0. 



if our first Mellin moments numerically fulfill bounds analogous to those derived for TMD 
PDFs in Ref. [BBHM00]. Consider the bounds on g\x given by eq. (2.57). In Fig. 4.16, we 
plot the function 

a,(i*ld -^(A (i)B >i)) 2 -(^r(fci)) 2 - 

with respect to |fej_| on the coarse-m020 ensemble for u and (f-quarks, and find that it is 
positive. Thus, on a numerical level, we find that the first Mellin moment <7^ sW (fc^_) complies 
with a similar bound as the corresponding ^-dependent TMD PDF. 

In the following, we want to give our Mellin moments an interpretation in terms of densities 
of quarks (minus antiquarks). In particular, we consider 

= / 1 (1) (fci), (4.86) 
PTL(k ± ; S ± , A) = 1 $[7 + KD (fc±; P, S) + * $[7 + 7 5 ](D (fc±; P, S )\ s=s± 

= lfi 1 \kl) + l k ^ 9 V(ki). (4.87) 
Here pTL(k±; S±, A) is defined in analogy to eq. (2.56). 

The unpolarized density pjju is the density of quarks (minus antiquarks) in a nucleon with 
respect to the intrinsic transverse momentum of the quarks, averaged over the nucleon spin. 
We plot this density in the two-dimensional transverse momentum plane in Fig. 4.17a. The 
distribution is axially symmetric. 

Particularly interesting is prL, the transverse momentum dependent density of quarks (minus 
antiquarks) with definite helicity in a transversely polarized nucleon, compare section 2.7.3. 



m N 



(4.85) 
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Figure 4.17: Quark density plots in the transverse momentum plane at m n = 500 MeV ob- 
tained from the Gaussian fits to the amplitudes A2(£ 2 , 0) and Aj(l 2 , 0) as depicted in Fig. 4.12. 
We highlight the origin with a white cross. 

(a) Density of up quarks (minus antiquarks) Pij(j u 

(b) Density p^f[ u (k±; S±, A) of up quarks (minus antiquarks) with positive helicity A = 1 (i.e., 
with spin pointing in ^-direction) in a nucleon polarized in transverse x-direction S± = (1,0). 
The upper left inset symbolizes the nucleon with its spin vector pointing to the right, and the 
probed quark with its spin pointing towards the reader. 

(c) Same as in (b) but for down quarks. 

(c) Same as in (b) for up quarks minus down quarks. The deformation appears amplified. 
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We plot this distribution for up quarks in Fig. 4.17b. For an intuitive interpretation, imagine 
a proton is moving towards us, the observer. The proton is polarized in transverse direction, 
with the spin vector pointing in positive x-direction, S± = (1,0). We now measure the 
probability distribution of u-quarks with spins pointing towards us (more precisely, with 
positive helicity A = 1). The density of quarks in the transverse momentum plane thus 
obtained appears deformed - it is no longer axially symmetric. The peak of the density is 
shifted to the right, to a positive value of fcj_i- The reason for this deformation is that gyr 
is non-zero. For the corresponding <i-quark distribution plotted in Fig. (4.17c), the peak is 
shifted to the opposite direction, a result of the sign change in Af or rather g^, see Figs. 
4.12 and 4.15. The deformation appears amplified if we plot pTL.(k±; S±, A) for u — d quarks 
as in Fig. 4.17d. 

Within the formalism of Ref. [Mil07], the deformations we observe in ptl are evidence for the 
"non-spherical shape of the nucleon" . The densities plotted in Fig. 4.17 are reminiscent of the 
spin densities in the spatial transverse plane presented in Ref. [G + 07] . In contrast to the study 
at hand, these densities have been obtained from GPDs, which have an interpretation in terms 
of impact parameter dependent distributions [DH05] . Concerning the comparison to impact 
parameter fa^-dependent densities, the distribution g^(k\) is somewhat peculiar. It describes 
a correlation between transverse nucleon spin, quark helicity and transverse momentum of the 
form XS±-k±. A corresponding impact parameter dependent distribution with a correlation 
of the form XS±-b± does not exist [DH05]. Another peculiarity of g^{k\) and its associated 
density ptl is the shift of quark density in the direction of the transverse nucleon spin. In 
contrast, the densities plotted in Ref. [G + 07] feature an asymmetry perpendicular to the spin 
vector, due to a correlation in the form of a vector product S± x b±. 

To give a quantitative description of our findings, let us calculate some fcj_-moments: 

2 \ _ / d2k ± k ± Puu{k ± ) 



(hi 
(fej 

9A 
9V 



Puu 



f d 2 k ± puu{k ± ) 
= / d 2 k± k ± pTL(k±;S±,X) 
} p™- fd 2 k ± p TL (k ± ;S±,X) ' 

fd*k ± &-r + ^)(kr,P,S)\ A=lt Si 



(4.88) 



TMD Jd 2 k± $[7 + ](l)(fc_ L ;P,S) 

Here gA/dv is a quantity that can be obtained without any reference to transverse momentum 
dependence. The nucleon vector coupling constant gy is given by the number of valence 
quarks (i.e., 2 in the u— and 1 in the d— or (u— cQ-channel). The nucleon axial vector coupling 
constant g^ for u—d quarks has been determined experimentally with high accuracy from 
neutron /3-decay: [gA/ 9v]^d = 1-2695(29) [A + 08]. Using the Gaussian parametrization of 
the amplitudes, we find 



Puu 

9_a] _C6 
.9v\tmd C 2 



171NC7 
C2 



(4.89) 
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Again, we remind the reader that the integrals in eq. (4.88) only exist if the distributions 
decay sufficiently fast with k±. With the Gaussian parametrization, this is guaranteed, but 
the perturbative high-|fcj_| behavior is not reproduced correctly. 

Numerical results for the fcj_-moments are included in Tables 4.3, 4.4a and 4.4b. Our results 
renormalized with the static quark potential and C ren = yield (k 2 L ) Puu « 0.4 GeV, which is of 
the same order of magnitude as the experimental estimate 0.5 GeV of Ref. [A + 05] mentioned 
in section 2.7.1. However, a serious quantitative comparison of our results for the RMS 
transverse momentum to phenomenological estimates is not justified at this stage, primarily 
for the following reasons: Firstly, in contrast to [<7a/<7v]tmd and {k±) PTL , the RMS transverse 
momentum (k\) puu is very sensitive to 5rh, and thus to the choice of the renormalization 
condition. Moreover, we have employed a simplified contour for the Wilson line and pion 
masses much larger than the physical ones. 

The numbers in Table 4.4b for (fcj_) pTL are in line with what we saw in Figs. 4.17b and 
4.17c: Helicity polarized quarks in a transversely polarized nucleon carry a non-zero average 
transverse momentum. This average transverse momentum shift differs in sign for u- and 
(f-quarks. Specifically, with renormalization from the static quark potential, we find a shift 
of (k±) PTL = (73 ± 5)ASj_MeV for up quarks, and a shift of about half the magnitude, 
(k±) PTL = (—31 ± 5)XS± MeV, for down quarks. 

For u — d quarks, we find \qaI Sv]tmd = 1.21 ±4 on the coarse-m020 ensemble with renormal- 
ization from the static quark potential, see Table 4.4a. We find agreement within errors with 
the value qaIsv = 1-173(29) determined on the same ensemble in Ref. [E + 06]. Notice that 
our result is not so far from the experimental value. A quantitative comparison to experiment 
can of course only be made with an extrapolation to the physical pion mass, as performed, 
e.g., in the above reference [E + 06]. We rate the successful determination of gA/dv as an 
important crosscheck of our methods. 



4.6.6 Ratios of Amplitudes 

Observables which can be obtained from ratios of amplitudes Ai are particularly attractive, 
because renormalization factors cancel entirely: 

Mt 2 ,Q) = ^r cn (i 2 ,Q) u 90) 

A 2 (£ 2 ,0) i™(^ 2 ,0) ' 

The unrenormalized amplitudes can be directly obtained from the unrenormalized ratios, 
provided I 2 is in the range where discretization artefacts are small enough. We plot ratios of 
amplitudes in Fig. 4.18. 

What can we read off from these amplitudes? The formalism in appendix B.6 shows how 
fc_i_-moments can be expressed in terms of the amplitudes A(£ 2 ,0) at I 2 = 0. Concerning the 
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(a) 



ens. 


quarks 


—5m 


x 2 


Cq 


<7fi (fm) 


9a/ 9v 


m020 


u — d 


1553(47) 


0.2 


-1 211(36) 


1 173(33) 


1 21lC36)d4) 


11 


u 


» 


0.5 


-0.927(33) 


1.202(39) 


0.458(16)(06) 


11 


d 


11 


0.9 


0.287(18) 


1.066(52) 


-0.281(18)(03) 


m020 


u — d 





1.9 


-1.124(34) 


0.7713(72) 


1.124(34) 


55 


u 


11 


1.8 


-0.860(31) 


0.7789(94) 


0.425(15) 


)) 


d 


11 


1.0 


0.269(17) 


0.734(19) 


-0.262(17) 


m030 


u — d 





4.3 


-1.149(20) 


0.7638(40) 


1.149(20) 


m050 


u — d 





9.6 


-1.123(13) 


0.7541(34) 


1.123(13) 



(b) 



ens. 


quarks 


—5rh 


x 2 


C7 


<r 7 (fm) 


(k ± ) PTL (GeV x XS ± ) 


m020 


u — d 


0.1553(47) 


0.4 


-0.1087(50) 


1.096(31) 


0.1793(83)(30)(63) 


)) 


u 


11 


0.6 


-0.0895(53) 


1.094(37) 


0.0730(42)(13)(26) 


11 


d 


11 


1.1 


0.0190(29) 


1.103(97) 


-0.0307(48)(06)(11) 


m020 


u — d 





0.7 


-0.0970(46) 


0.7586(95) 


0.1600(76) 


)? 


u 


11 


0.8 


-0.0796(47) 


0.760(12) 


0.0649(38) 


)) 


d 


11 


1.1 


0.0169(27) 


0.765(36) 


-0.0271(44) 



Table 4.4: Numbers obtained from the Gaussian fits for Aq and Aj, determined on the HYP 
smeared coarse ensembles. The value for 5m in the upper sections has been determined using 
the static quark potential, see section 4.4.7 in the text. The renormalization constant is 
taken from the analysis of A2 jU -d- Due to correlations, the value \ 2 (P er degree of freedom) 
specified here gives only a crude orientation about the fit quality, see section 3.4.3. The 
meaning of the errors is, in sequence, (1) statistical error, determined from 1000 bootstrap 
samples, (2) crude estimate of systematic errors from the smeared-unsmeared comparison of 
section 4.4.6, and (3) uncertainty from the determination of the lattice spacing. 
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Figure 4.18: Ratios of amplitudes Ai{£ 2 ,£ ■ P = 0) to A 2 (£ 2 ,£ ■ P = 0) for straight Wilson 
lines calculated on the HYP smeared coarse-m020 ensemble. These quantities need no renor- 
malization. The shaded regions in the range V^P < 0.25 fm and for V^I 2 > L/2 give a 
conservative estimate where we expect lattice cutoff effects and finite volume effects. 
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ens. 


quarks 


9A/9V 


[(k±) PTL ] reg (GeV x XS ± ) 


m020 


u — d 


1.216(48) 


0.162(19)(06) 


5) 


u 


0.468(21) 


0.072(10)(03) 


11 


d 


-0.256(23) 


-0.017(12)(01) 


m030 


u — d 


1.199(30) 




m050 


u — d 


1.193(22) 





Table 4.5: Numbers obtained from the renormalization scale and scheme independent ratios 
A1/A2, determined on the HYP smeared coarse ensembles. The first error is statistical and has 
been determined using 1000 bootstrap samples, the second error comes from the determination 
of the lattice spacing a. Unspecified systematic uncertainties also arise from the polynomial 
extrapolation of the ratios of amplitudes down to £ 2 = 0. 



observables of the previous section, we readily obtain 

_ d 



Puu 
9A 



d(-P) 4(^ 2 , 



\e=o 



4(0,0) 



9v 
(k±) 



4,(0,0) 



PTL 



TMD 4(0,0) ' 

= _ xs m N A 7 (0,0) 

reg 1 i 2 (0,0) 



(4.91) 
(4.92) 
(4.93) 



In eqns. (4.91) and (4.93), the subscript "reg" means that we have omitted certain terms that 
vanish if we make the additional assumption that derivatives of Ai{£ 2 , 0) are finite at £ 2 = 0. 
The omitted terms are higher derivatives of the amplitudes multiplied by factors of £± or £ 2 , 
compare the right hand side of eq. (B.23) in the appendix. Incidentally, the assumption of 
finite derivatives at £ 2 = is true for Gaussian amplitudes, so that we recover eq. (4.89) upon 
substitution of the Gaussian parametrizations. 

The solid curves and error bands in Fig. 4.18 have been obtained from fits to polynomials of 
second order, where only the data for 0.25 fm < \J — £ 2 < L/2 has been used as input, and 
where the fit weights have been chosen according to the statistical errors of the individual data 
points. Let us make an extrapolation of the ratios Ai(£ 2 , 0)/A2(£ 2 , 0) down to £ 2 = based 
on the polynomial fits. This way, we obtain numbers for [9a/ 9v]tmt> and [(k±) pTi ] reg , which 
we list in Table 4.5. The numbers qualitatively agree with the results from the Gaussian fits. 
Small differences are expected, because the ratio method and the Gaussian parametrization 
represent two different ways to extrapolate down to £ 2 = 0. In that sense, these differences 
give a crude estimate for a systematic uncertainty. The numbers for [(m/svJtmd in the (u—d)- 
channel agree within their (relatively large) statistical errors with the results for quoted 
in Ref. [E + 06], although we note that they are systematically higher. 



We now understand why gAj 9v and (fej 



'PTL 



turned out to be insensitive to 5fn in the previous 



chapter: They can be obtained from ratios of amplitudes. It will be interesting to study the 
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Figure 4.19: Linear chiral extrapolation of the unrenormalized RMS width of the transverse 
momentum distribution (k 2 t ) obtained from the coarse ensembles. The error bars and the 
error band correspond to the statistical errors. We indicate the location of the physical pion 
mass by a dashed vertical line. 



small-£ 2 -behavior of the amplitudes and its effect on the quantity {k±) PTL . 



4.6.7 Chiral Extrapolation 

To make physical predictions, an extrapolation to physical quark masses is needed. For 

some quantities, such as qa, sophisticated extrapolation formulae have been worked out in 

chiral effective field theory. Using a simple linear extrapolation instead introduces large 

unquantifiable systematic errors but may serve as a first guess. Let us look at quantities 

specifically related to TMD PDFs. With the data presently available from our analyses, 

we can only investigate the quark mass dependence of the RMS width of the transverse 

momentum distribution (k 2 t ) as we obtain it from the Gaussian fit. Up to now, we 

x ±/ Puu n 
have not worked out the renormalization constant 5rh for the coarse-m030 and coarse-m050 

ensembles, but it is plausible that they will not depend strongly on the sea quark masses. 
Let us take the values for u—d with 5rh = from Table 4.3 and extrapolate them 

linearly down to the physical pion mass. We obtain (fej.) 1 ^ = 0.5558(94) sta t GeV (with a x 2 
per degree of freedom of 0.7). The extrapolation is illustrated in Fig. 4.19. It turns out that 
the extrapolated value hardly differs from our input at large m n . For the high pion masses in 
our study, the RMS transverse momentum shows practically no sensitivity to the light quark 
masses. However, we should be aware of the fact that a much more pronounced quark mass 
dependence may set in at lower pion masses 8 . 



3 as is expected in, e.g., the extrapolation of g A [HPW03, K+06, PMHW07, E+06] and of (x) [DGH07]. 
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10 12 14 




(fin) 



Figure 4.20: Illustration of our procedure to obtain cuts at constant £ 2 . The black dots are 
points in the (£ 2 ,£-P) -plane where lattice data is available. The vertical lines are selected 
values of £ 2 for which we plot A2(£ 2 ,£-P) in Fig. 4.21. The colored dots correspond to the 
lattice data points shown in these figures. To obtain intermediate values, we have interpolated 
linearly on the (£ 2 , £-P)-plane, using the triangulation shown as a grey mesh. 



4.7 Dependence on the Longitudinal Momentum Fraction x 



In the previous section we have studied lattice data for £-P = 0. Let us now explore the 
(^•P)-dependence of A2(£ 2 ,£-P). The first integral in the last line of eq. (4.17) shows that it 
is related to the dependence of f\ at (x,k±) on the longitudinal momentum fraction x via a 
Fourier transformation. 

We gave a three-dimensional overview of the available data for this amplitude in Fig. 4.4. Let 
us now work on the coarse-m030 ensemble, where we have rather good statistics. In order to 
be able to plot A 2 (£ 2 ,£-P) as a function of £ • P at fixed values of £ 2 , we interpolate linearly 
between our data points as illustrated in Fig. 4.20. This enables us to take a closer look 
at the (^-P)-dependence in Fig. 4.21. First of all, the symmetry of the plots clearly confirms 
the relation [Ai(£ 2 , £-P)\* = A{(£ 2 , —£-P) we found in section 4.1.1. Moreover, we recognize 
a general suppression of the amplitude with increasing £ 2 . Apart from this suppression, and 
apart from finer structures which might be statistical fluctuations or discretization errors 
from the steplike linkpaths, the (£-P)-behavior follows a general trend independent of £ 2 : 
In the real part, we notice a slight curvature downward. The imaginary part shows strong 
(£-P)-dependence in the form of an almost linear rise. 
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Figure 4.21: Cuts through the unrenormalized amplitude A2(£ 2 ,£-P) for a selection of fixed 
values y/—£ 2 = \£\ on the smeared coarse-m030 ensemble for u — d quarks. The error band 
between the data points has been obtained from linear interpolation as illustrated in Fig. 4.20. 
(a) Real part, (b) Imaginary part. For the sake of clarity, we have added offsets in the 
ordinate in steps of 0.1. The dashed lines indicate the respective zero lines. 
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4.7.1 A Normalized Amplitude 

The observations above have inspired us to study the normalized amplitude 

Ar m (i 2 ,e-P) = Me2j ' p) = j|!^Vjg) (4 .94) 

2 V ' A 2 (t 2 ,0) Af licn (t 2 ,0) V ' 

We point out that the quantity defined above needs no renormalization, just as the ratios of 
amplitudes in eq. (4.90). At present, we obtain A^ orm (i 2 , t-P) from 

n I p /->lat\ 

A%°™(t 2 ,t-P) = ^(P) ■ (495) 

Here i? 73 (P, C] at ) is the (unrenormalized) ratio defined in eq. (3.58). The function I(i 2 ) is 
a smooth interpolation of Ry i (P, C^ at ) at £-P = 0, where we have chosen straight link paths 
Ci at on the lattice axes. In the future, we should use the ratio 



AT vm U 2 ,i-P) = Ti{ ,V - (4.96) 



The expected advantage of the latter is an optimal cancellation of statistical fluctuations and 
residual discretization errors from the step-like link paths. Up to now, data for i? 7a (0,C£ at ) 
has not been calculated for the full set of link paths, so we must content ourselves with the 
interpolation method. 

We plot A2° Tra (£ 2 , i-P) as a function of t-P for fixed £ 2 in Fig. 4.22. It appears as though 
the normalized amplitude is largely £ 2 -independent. In the real part, we observe a slight 
curvature; the amplitude bends down with increasing \£-P\. A much cleaner picture emerges 
for the imaginary part. Here the amplitude basically describes a straight line through the 
origin, with a slope that appears to be ^-independent. 

The values for i-P we can access on the lattice are integer multiples of 2n/L. In Fig. 4.23 we 
plot A2° rm (£ 2 , i-P) for each of these values as a function of i 2 . We do not see any significant 
dependence of A2° rm (i 2 , i-P) on i 2 . In conclusion, within the range of available lattice data 
and within our level of precision it is appropriate to write 

I| orm (£ 2 , i-P) w A% oria (£-P) (4.97) 

and, making use of definition eq. (4.94), 

A 2 (i 2 , i-P) » A^ii-P) A 2 (i 2 , 0) . (4.98) 

In other words, the amplitude A 2 (i 2 ,i-P) approximately factorizes into an £ 2 - and an (^-in- 
dependent part. 



4.7.2 Factorization Hypothesis 



What might be the implication of this observation in terms of TMD PDFs? Suppose eq. (4.98) 
were true in the entire domain of £ 2 and t-P. Inserting this hypothetical factorization into 
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Figure 4.22: Cuts through the normalized amplitude 
values \/—£ 2 



l (£ 2 ,£-P) for a selection of fixed 



\£\ on the smeared coarse-m030 ensemble for u — d quarks. The error band 
between the data points has been obtained from linear interpolation as illustrated in Fig. 4.20. 
For the sake of clarity, we have added offsets in the ordinate in steps of 0.5. 

(a) Real part. The dashed lines are plotted at an ordinate of 1, shifted by the respective 
offset. 

(b) Imaginary part. The dashed lines indicate the respective zero lines. 
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Figure 4.23: Testing ^-independence on the smeared coarse-m030 ensemble for u — d quarks: 
We plot our data for Af Tra {l 2 ,£-P) with offsets (£-P)L/2ir in the ordinate. To guide the eye, 
we display the average of A2° rm (£ 2 ,£-P) at constant l-P as horizontal dashed lines. In the 
averaging procedure for the horizontal lines, we also include data points with l-P of opposite 
sign, taking into account the property [i^ orm (£ 2 , l-P)]* = Af Tra (l 2 , -l-P). The magnitude 
of the amplitude relative to the offset can be read off from the inclined dashed lines on the 
left. The regions with a gray background are more likely to be affected by significant lattice 
artefacts. 
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eq. (4.17), we obtain 

7 27T A 2 (0,0) Jo 2 ( 27r ) 

V v ' ' v ' 

/i(x)/i 2 (0,0) tf\k x ) 

(4.99) 

The two integrations can now be carried out independently, so that the factorization assump- 
tion of A2(£ 2 ,£-P) translates into a factorization of the quark density 

a Mx,k±) = ^fP(k±) n (4.100) 

into an x- and a fcj_-dependent part. The latter is just the first Mellin moment discussed 
in section 4.6. The x-dependent part is identified as the regular PDF f\(x) divided by the 
number of valence quarks M = ^2(0, 0). 

At this point, some serious words of caution are in order. The observations we make with lim- 
ited statistical significance and within the limited range of £ 2 and £-P on the lattice certainly 
do not justify a claim that A2(£ 2 ,£-P) factorizes strictly in the entire domain. Moreover, the 
amplitude A2° rm (£ 2 , £-P) needs no renormalization, while, in contrast, fi(x) is a renormal- 
ized quantity. Obviously, somewhere in our calculation eq. (4.99) we have ignored necessary 
renormalization steps. As discussed in section section 4.6.3, we expect that in the contin- 
uum the amplitude A2(£ 2 ,£-P) diverges for £ 2 — » 0. Thus we may reason that the expression 
A2(0,£-P)/A2(0,0) is only well defined within a suitable regularization prescription, which 
then carries over to f\{x). 

At the present stage, eq. (4.99) shows us that there is a qualitative analogy between the 
observed factorization of ^4 2 (^ 2 , £-P) and a factorization assumption of the form of eq. (4.100). 
The latter has been frequently used as a working hypothesis in phenomenological applications, 
in particular in combination with the Gaussian parametrization of /{ (fcj_), see, e.g., the 
analysis of Ref. [A + 05] discussed in section 2.7.1. Our lattice data provides some evidence that 
the factorization assumption is justified as an approximation, at least in a certain kinematical 
region, with the caveat that we have used straight Wilson lines in our correlators. 

In Fig. 4.24, we plot our lattice data for A2° rm (£ 2 , £-P) in stripes of constant £ P, omitting 
data with \/—£ 2 < 0.25 fm to be on the safe side concerning lattice cutoff effects. Small offsets 
on the abscissa enable us to verify again that there is no strong £ 2 -dependence at a given £-P. 
The error bands superimposed on the data have been obtained with polynomial fits which 
correspond to a parametrization of A 2 lorm as 

A™™(£ 2 , £-P) = 1 + i(£-P) b 2 - l(£-P) 2 63 - k£-P) 3 b 4 + ^.(£-P) 4 b 5 . (4.101) 

2 6 24 

Again, the fit weights have been chosen according to the statistical errors of the individual 

data points. The values we obtain from our fits in the u — d channel (where jV = 1) are 

b 2 = 0.242(19) , 63 = 0.226(49) , 64 = 0.019(18) , 65 = 0.157(60) . (4.102) 
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From the factorization assumption eq. (4.99) follows 

dx x n fi(x) = n J *£p. ir-a-p) e-^ 



M(-i) n ^J^A° Tm (e-P) =Mb n+l . (4.103) 



d(£-P) n 



■p=o 



Thus, if we could take eq. (4.99) literally, the b n would be Mellin moments of f\{x). An 
analysis of x-moments based on renormalized local operators on the same ensemble [H + 08] 
finds J dx x fi, u -d( x ) = 0.226(4), which is quite similar to our however, our 63 differs 
strongly from their result f dx x 2 f\ )U -d{x) = 0.074(5). 

4.7.3 Qualitative Comparison to PDFs from Phenomenology 

Taking eq. (4.99) literally once more, we can invert the Fourier transformation and get 



" A^(e-P) = AT 1 / dx eW p > f 1>q (x) 



1 

-1 



= AT- 1 f 1 dxe^ p > i f M X K : ", (4.104) 

y_! 1 -h,q\- X ) ■ X < 

where we have made use of the support properties of fi(x). Let us take the phenomenologically 
determined CTEQ5M parametrization [L + 00] of fi tU -d(x) and f 1 u-d( x ) a ^ a scale Q 2 = (2 GeV) 2 
as input to the above equation. An unofficial release of the CTEQ5M distributions is available 
as a Mathematica input file. We have evaluated the above integral numerically, excluding a 
small x-interval [—0.00001,0.00001] to avoid numerical problems. The results are shown as 
red dashed curves in Fig. 4.24. At a qualitative level, there is an obvious similarity to the 
lattice data. 



4.7.4 Qualitative Comparison to the Diquark Model 

The model of section 2.8 provides us with an explicit expression for the quark distribution 
fi(x,k±) valid for x > 0. We can convert this into an amplitude A2 via the inverse Fourier 
transformation 

rl poo 

A 2 (£ 2 ,£-P)= dxeW p > d\k ± \ 2vr|fc ± | Mv^i 2 \k ± \) h(x, k\) (4.105) 
Jo Jo 

provided fi(x, k\) decays asymptotically fast enough with fcj_ (otherwise, some regularization 
procedure is required, e.g., via a cutoff). Let us take the parameter set from section 2.8 in its 
original form with a = 2, as in Ref. [JMR97]. In this case, we can evaluate the above integrals 
numerically without a cutoff in From the resulting amplitude A2(£ 2 ,£-P), we calculate 

A2° rm (£ 2 ,£-P) according to its definition eq. (4.94). We plot it in Fig. 4.24 with respect to 
l-P for two different values of £ 2 , namely \/—£ 2 = and \J — £ 2 = 1 fm. 
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First of all, we find qualitative agreement between the model and lattice data. What is 
perhaps even more astonishing is that the model amplitude JlJ} ™ 1 ^ 2 , i-P) is very similar for 
\/—£ 2 = and y/—£ 2 = 1 fm. We conclude, that the factorization assumption on the level of 
amplitudes A2(£ 2 ,£-P) is also a good approximation for the scalar diquark model within the 
region of (£ 2 ,£-P) accessible to us on the lattice. 



4.7 Dependence on the Longitudinal Momentum Fraction x 
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Figure 4.24: Results for A^ orui {i 2 , £-P) obtained on the smeared coarse-m030 ensemble for 
u — d quarks. The statistical error bars of the lattice data are shown as colored rectangles; 
data with larger statistical errors are drawn with lighter colors. In each vertical stripe we show 
lattice data for one value of £-P; data points obtained for larger values of yj—l 2 are drawn 
further to the right inside the stripe. We have omitted data points with y — £ 2 < 0.25 fm 
because of possible lattice cutoff effects. The curves labeld "CTEQ" and "diquark model" are 
explained in the text and are only intended to provide a qualitative comparison. 
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In this chapter, we would like to present some preliminary studies and ideas that may help us 
in the future. In particular, we will address the question whether we can go beyond straight 
Wilson lines in our operator, with the aim to calculate "realistic" TMD PDFs, as they occur 
in the description of scattering experiments such as SIDIS. 

5.1 Auxiliary Fields 

The operator q{£) r op K[Cg] g(0) is non-local. This complicates the probability interpretation 
and renormalization of our quark-quark correlator. Therefore, it can be useful to give the 
gauge link an alternative interpretation as a propagator of some auxiliary field, which produces 
the gauge link when the auxiliary field is integrated out. We draw attention to the fact that 
such an auxiliary field formalism ("z-field formalism") has already been used in the literature 
to derive the renormalization properties we make use of in section 4.4.1. Moreover, the idea 
of integrating out a fermion is somehow linked to the motivation for the introduction of the 
Wilson line in the first place: As we discussed section 2.4.2, the longitudinal Wilson line 
effectively describes interactions with a fast propagating parton. In the following, we show 
why auxiliary field techniques (and effective theories) can become interesting for us. 

5.1.1 Heavy Particles 

That heavy auxiliary fields can be advantageous in the context of hadron structure calculations 
on the lattice has already been put forward in Ref. [DL06] . The heavy quark action is derived 
from the QCD fermion Lagrange density for a single heavy quark flavor: 

C F [¥, A] = ^{x){ilp - m h )^{x) . (5.1) 

Here is the mass of the heavy quark field Due to its large inertia, the motion 

of the heavy quark is characterized by small fluctuations around the classical path. Let us 
decompose its momentum according to p = uiqv + k, where v is a fixed velocity four vector 
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(v-v = 1) and where tuq ~ m^. The choice of uiq and U are arbitrary to a certain extent, as 
long as the remaining dynamical momentum k is small compared to tuq. One now decomposes 
^f(x) such that 

1 + <& 1 — 4 

— - — ^(x) = exp(—imQ v-x) h(x) , — - — ^(x) = exp(— irriQ v-x) H[x) . (5-2) 

with projections (1 db f)/2. The phase factor absorbs oscillations due to the motion with 
velocity v. Making use of the classical equations of motion for the field H(x), we arrive at 
the HQET Lagrangian [Geo90, Neu94] 

£hqet[/*, h, A] = h{x) (iv-D - - m^j h(x) + 0{l/m%) . (5.3) 

Here itir = rrih — rriQ is the residual mass [FNL92] and Dj, = — {v-D)v^. Choosing v in 
temporal direction and going to Euclidean space, we get 

£% R0 [h, h, A] = h(x) (-D- 4 - ^ + m^j h(x) , (5.4) 

where we have omitted a spin dependent relativistic effect which is part of Ipx^Ti an d where 
D refers to the three spatial components of D. In the static limit mg — ► oo, the second 
derivative term vanishes. A simple discretization of the static action on a lattice L of infinite 
volume then reads 

-Sstatt^ h ,U] = Y^ H x ) [K x ) - U{x, x + l)h(x + 4)] + h(x)m R h{x) . (5.5) 
A solution for the propagator of a static quark is thus 



h(x)K{y) = (1 + mfl) w - Xi_1 



U(x, x + 4)U(x + 4, x + 2 4) • • • U(y - 4, y) 

1 





x j — Vj: 2/4 > x 4 

x = y ■ (5.6) 

otherwise 



The quark can only propagate in one direction in Euclidean time. The propagator becomes 
a straight gauge link in 4-direction, once the quark fields are integrated out [E1C88]. 1 We 
have already alluded to this formalism in section (4.4.7.1), when we interpreted the temporal 
line segments of Wilson loops as propagators of static quarks. Note that (1 + rhR) Vi ~ Xl ~ 
exjp(mji(yi — xg)). This factor resembles the Wilson line self energy exp(— 6my/— £•£). 

For our quark-quark correlator with the straight Wilson line, we would like to find an action 
that produces the straight link between any two points separated by a spatial vector I. Ideally, 
this action will be independent of i. In pursuit of such an action, let us study eq. (5.4) again. 
This time, we also take the term D 2 /2mg into account. We remark that the whole expression 



x We set the fermion determinant of the heavy quark action to one ("quenched approximation"), i.e., we 
neglect virtual heavy quark loops. 
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eq. (5.4) is just the leading order Lagrangian in the NRQCD formalism [CL86, BBL95], which 
has a different power counting scheme than HQET (compare, e.g., Ref. [LM97]). A simple 
discretization is 

1 3 

= ^^W,y)%) , (5.7) 

where 

K(x,y) = S x , y (l + m R + — J -$ x>y -iU(x,y) - — - ^ 5 x>y -p,U(x,y) . (5.8) 

^ c— 

9 =M(x,y) 
i — t~i 

The propagator G(z,x) = h{z)h{x) fulfills 



5 z , y = G(z, x)K(x, y) = P G(z, y) - G(z, y - 4)U(y - 4, y) 

I- G(z,x)M(x,y) . (5.9) 



2mr, 



Let us study the restriction of the above equation to z\ = y\. If we do not allow the particle 
to propagate backward in Euclidean time (analogous to eq. (5.6)), then G(z,y — 4) = 0. We 
can now solve for G (in matrix notation): 



2m Q J pj^ Q \ 2 P™Q J 

Now, consider the propagator of the particle from the origin to £, with = 0. Higher powers 
of M are suppressed by powers of 1/ttiq in the series above. The matrix M connects a given 
lattice site with the six neighboring lattice sites in the same time slice. The expansion in 
powers of M is called a hopping parameter expansion, see, e.g., Ref. [Rot97]. The leading 
non- vanishing contribution to the propagator has [Y| = Y%= l \^p\ factors of M: 

h(£)K{0) = p- 1 (2pm Q )-W \ ^M(x,z(W-i))...M(z( 2 )z«)M(zW0) + 0(m Q 1 ) 

Z (i),...,«(m-D 



lc lat 



C lat to I 



+ 0{m 1 ) \ . (5.11) 



If i is chosen to lie on one of the lattice axes, we immediately see that the leading order 
contribution is the shortest possible gauge path from to b. a straight link. This leading 
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contribution comes with a factor ~ tHq ; again this factor reminds us of the Wilson line 
self-energy exp(— <5m|^|). The next order introduces a dent in the path, see Fig. 5.1a. The 
dent is suppressed by 0(toq 2 ), but there are many possible locations along the path where 
the dent can be inserted, and the dent can have different lengths. If the vector £ is at an 
oblique angle compared to the lattice axes, many link paths contribute to the leading order, 
see illustration Fig. 5.1b. Link paths running close to the geometrically direct connection 
dominate combinatorially, as illustrated in Figs. 5.1c and 5. Id. Comparing the two figures, 
we note that the resemblance to a direct connection intensifies on a finer lattice. However, the 
continuum limit a — > oo cannot be taken within our picture, because the hopping parameter 
expansion is only applicable (probably as an asymptotic series) on a not too fine lattice, where 
arriQ S> 1. Nevertheless, we conclude that at leading order in itIq 1 for not too small \£\ the 
propagator effectively approximates a straight Wilson line from to £? It seems that for a 
given rriQ, this Wilson line has a certain "thickness", in the sense that gauge links with small 
dents deviating from the direct path contribute (no matter how fine we choose the lattice). 
The larger we choose mn, the "thinner" the line will get. In the continuum formalism, paths 
with tiny wiggles can be absorbed into the renormalization constant 5m of a straight Wilson 
line, see Ref. [CD81], so we see another correlation between tuq and 5m. In summary, we 
rewrite eq. (5.11) symbolically as 

" /iSo) « U thick [£,0] m Q W ", (5.12) 

where U thlck [£, 0] represents the combination of Wilson lines from to £, and where the factor 
tUq reminds us that there is a suppression factor that exponentiates with the distance \£\. 

How can we make use of these observations? With respect to the action S[q, q, U]+S^^o[h, h, U], 
we can consider the correlator 

$ h (kx,P,s) = J ^0^e- ik -^(N(p,s)\ m^Ht) Mo)rVo) \N(P,S))\ £±=0 

(5.13) 

(in analogy to the first Mellin moment in eq. (2.13)). Here the Dirac structure r op of our 

original correlator eq. (2.5) is encoded in T h . More research is needed to understand the role 

of T h , and the role of spin in general in this context 3 . We now integrate out the heavy degrees 

of freedom. This replaces h(£)h(0) in the expression above by the heavy quark propagator 
i — i~i 

h(£)h(0). As we have demonstrated, this propagator translates into a Wilson line, at least in 
an approximate sense. Within the symbolical notation of eq. (5.12), we thus obtain 

" * h (k ± ,P,S) « J e-^ \ (N(P,S)\ q(£) T^T h 

xU thick [£,0] m Q W q(0) \N(P,S)) \ e+=E . =0 "• (5.14) 

2 Note that our "leading contribution" might not be the leading one if we do not truncate our heavy quark 
action in the beginning. 

3 Note that the heavy quark fields h(x) involve a spin projection, see eq. (5.2). This restricts valid choices 
ofT h . 



5.1 Auxiliary Fields 



129 



(a) (b) 




Figure 5.1: Hopping parameter expansion for a spatial propagator in our heavy particle action, 
(a) Leading link path and example of a subleading link path for propagation along a spatial 
lattice axis, (b) Examples of leading link paths for propagation at an oblique angle, (c), (d) 
Superposition of all leading link paths for propagation at an oblique angle. The line thickness 
of a each link variable is proportional to the number of times the respective link variable 
occurs in the whole set of link paths. 
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Similar as in section 2.5.2.1, but without the difficulties with the gauge link, we can rewrite 
eq. (5.13) as 



for a complete set of states \n). This is the probability that a fictitious interaction converts a 
quark into an /i-particle such that the final state |n) carries transverse momentum — k±. 

5.1.2 How Auxiliary Fields may help us 

The preceding section shows that we can build correlators that look quite similar to our 
original one in eq. (2.5) once the auxiliary fields are integrated out. However, heavy particles 
are not the only type of auxiliary particles that might be interesting to us. Soft collinear 
effective theory (SCET), for example, describes very fast rather than very heavy quarks, 
and can generate the Wilson lines out to infinity that appear in SIDIS factorization. In 
general, the auxiliary field and action need not have a direct interpretation in the context 
of the actual physical process we want to study (e.g., SIDIS). The primary goal is to design 
well-defined observables that provide us with information about the transverse momentum 
distributions of quarks inside the nucleon. As demonstrated above, such an observable can 
have a probability interpretation. Moreover, we may learn more about the meaning of the 
renormalization condition. Last but not least, we might be able to construct improved lattice 
operators with reduced discretization artefacts, profiting, e.g., from existing experience with 
heavy quarks on the lattice. 



The quark-quark correlator defining TMD PDFs in the context of scattering experiments 
contains a Wilson line with sections running close to the ra_ direction, see chapter 2. Can we 
construct such correlators on the lattice? 

5.2.1 Accessible Link Directions on the Lattice 

Let us first visualize the situation. Suppose the operator we want to use to probe the nucleon 
contains a straight Wilson line whose direction is given by the Minkowski vector v. We want to 
implement this Wilson line directly on the lattice. In a frame of reference at rest relative to the 
lattice, v cannot have a Minkowski time component, v? t 0. (Neither do we consider Wilson 
lines with an extent in Euclidean time direction, because they complicate the interpretation 
in terms of the transfer matrix formalism.) For the purpose of our discussion, let us restrict 
ourselves to a vector of the form wi at = w sin(#)ni + w cos(#)n3, and consider a lattice nucleon 



$ h (k ± ,P,S) = -^J2 I < n l ^(0)rN(0) \N(P,S)) 2 6W(p n± + k ± ) >0 (5.15) 



n 
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Figure 5.2: Boosting the frame at rest relative to the lattice to the large momentum frame. 
In this illustration, the frame of reference is fixed by the requirement that P + = c l\f 7 lm^. 
The vector labelled "P" depicts the nucleon momentum. The other vectors represent v for 
= 0°, 30°, 60°, 90°, 120°, 150 °, and 180°. The dashed hyperbola e connect poin ts of equal 
proper length sj — (x ) 2 + (x 3 ) 2 = w and of equal invariant mass \J (x ) 2 — (x 3 ) 2 = mjv- (a) 
Nucleon at rest on the lattice, P^ t = 0. (b) Nucleon moving with P^ at = /2 with respect 
to the lattice. 



The quantities we want to calculate on the lattice have an interpretation in a frame of reference 
where the momentum of the nucleon is very large, so let us boost in 3-direction to a fixed large 
nucleon momentum P + . Figure 5.2 shows what happens to the vector v under this boost. For 
illustration purposes, we have chosen P + = 2\f2m^ (i.e., a Lorentz factor 7 = 2 with respect 
to the nucleon rest frame). In Fig. 5.2a, we start out with a nucleon at rest on the lattice, 
P^L = 0, while in Fig. 5.2b, we set Pj?+ = 771^/2, which is roughly comparable to our situation 
with the LHPC sequential propagators, where \P\ ~ 500 MeV. Under the boost, the nucleon 
momentum vector moves to the right along the dashed hyperbola at the top. Other vectors 
(except for purely transverse ones) are also pulled towards the h + axis. Of particular interest 
is the vector v for cos(#) = 1, whose arrow head touches the dashed hyperbola to the right. 
Obviously, by changing P^ t , we can move the arrow head up and down along the hyperbola. 
If P| 3 t is of the order of P + / \/2, the arrow is approximately horizontal, parallel to the B3 axis. 
If we were to increase Pj 3 t to a value of the order (P + ) 2 /mjv, the arrow head would move all 
the way down and appear close to the n_ axis. Of course, we can never reach Dtxn_ exactly. 

The direction v of the Wilson line relative to the nucleon momentum can be described by 
the Lorentz-invariant quantity £ = (P • v) 2 /\v 2 \ ~ (P + ) 2 /2 we already encountered in section 
2.4.3. In the limit of lightlike v, £ becomes infinite. On the lattice, £ is limited by the largest 
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(a) 



(b) 



TfO 



1 




Figure 5.3: (a) Staple shaped gauge link, (b) Wilson loop subtraction factor. 



attainable lattice nucleon momentum: < |Pi a t|. 



5.2.2 Staple-Shaped Wilson lines on the Lattice 

We have performed some preliminary studies on the lattice with staple shaped gauge links 
of the form Cg v = U[£,rjv + £, rjv, 0] as depicted in Fig. 5.3a. Again, we have implemented 
them directly as a product of link variables W lat [C] a * J . We have chosen a quark separation 

in 2-direction, I = £ 2 2\ and a staple vector v in 1-direction, v oc 1. In combination with 
the lattice nucleon momentum P = (2tt/L)(— 1, 0, 0), this choice enables us to work at non- 
zero v-P. Let us keep v normalized to v 2 = —1. With the replacement k — > imj^i, we 
deduce from eq. (7) of Ref. [GMS05] that for our particular choice of v and £, the ratio R r ™ 
should depend on the amplitudes A 2 , A\ 2 and B$ for r\ — > 00. At finite 77, let us denote the 
amplitudes by ai(£ 2 ,£-P, rjv-£, —rf, rjv-P) and bi(£ 2 ,£-P, 7]v-£, —rj 2 , rjv-P). In analogy to section 
4.1.1 and eq. (4.7), we obtain the ratio R™* 1 via structures A4. For our test case, the relevant 
.A/f-structures read 



M^(£, P, r,v) = — a 2 P» + 2i a 12 e^ af3 P a £p 7,7 s + 2im 2 N b 8 e^£ aV v p 7,7 s + 

Tfl j\r 



(5.16) 



In the limit of large rj, we expect the structures of the above equation to converge to finite 
values. The amplitudes we want to extract are then 

^2 (£ 2 ,£-P,j^,C\sgn(vP) ) = lim a 2 (£ 2 ,£-P, V v£, - V 2 , V vP) , (5.17) 

V \ V 'P\ J t]-^oo 

A 12 (l 2 ,£-P, T ^,C 1 ,sgn(v-P)) = lim ~a 12 {£ 2 ,£-P,r ] v£,-ri 2 ,i 1 vP) , (5.18) 

V \ V 'P\ J T)-^<X> 

Bs(£ 2 ,£-P,T^,C\sgn(vP)) = lmi(7]vP)b 8 (£ 2 ,£-P,riv£,-r] 2 ,7]vP). (5.19) 

\ \ V 'P\ J V^co 

Note that b$ must vanish at least as fast as 1/r] to keep i?g finite. Taking into account the 
transformation properties of the staple-like like link path eq. (2.48), the constraints on A4 
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from discrete symmetries eq. (4.4) must now be generalized to 



(t) : [Mr(i,P,v)f = 7 ° M jOV f^o(-£,P,v) 7 ° , 

0?) : Mr(£,P,v)= 7 ° M 7 o rj0 (I,P,v) 7 ° , 

CO = [Mr(l,P,v)]* =Y°C M ch s ri5c (-lP,-v) CV . 

From hermiticity (f ) follows for the amplitudes Ai (and analogously for the Bj) 

v£ 



\v-P\ 



The time reversal operation (T) implies 
v-l 



i-P, 



\v-P\ 



,( \sgmVP) 



A 2 

Al2 

Bx 



\v-P\ 

?,£-P,^,C\sgn(vP) 

?,£-P,^,C\sgn(vP) 
\v-P\ 



A 2 

-Al2 

-B» 



^ 2 ,-£-P,^,C 1 ,sgn(vP) 

,-£-P,^-,C\-sgn(v-P) 
\v-P\ 



(5.20) 
(5.21) 
(5.22) 

(5.23) 



■P,T^,C\sgn(vP)) . (5.24) 
\v-P\ J 



In combination with eq. (5.23), this confirms that A 2 is a T-even amplitude, A%(. . . , 1) = 
A 2 {...,-1), while Ai2 and B 8 are T-odd, Au(...,l) = -Au(. . . ,-1) and B&(...,1) = 
-B$(. . . , -1), as stated in Ref. [GMS05]. 4 

In our specific case, v-l = and £-P = 0, so according to eq. (5.23), all our amplitudes will be 
real valued. From our master formula eq. (3.61) and with the LHPC choice of spin projectors, 
we obtain the ratio 



R™(P,£,r)v) = 2a 2 



2imjv^2 
-EJPJ 



Pi an + m N r]vi 6 { 



(5.25) 



Our test calculation has been carried out on the 425 HYP smeared gauge configurations of 
the coarse-m050 ensemble. In Figure 5.4a, we plot the real part of the unrenormalized ratio 
Ry i (P, C^„) for our choices of £ and v. Clearly, the ratio goes to zero for large \i]v-P\. By 
now we know that this is the expected behavior: the gauge link on the lattice decays with 
increasing length. Partly, this decay can be attributed to a renormalization factor exp(2<5m|?y|). 
We realized that we can cancel all renormalization factors introduced by the gauge link if we 
divide by the square root of a vacuum expectation value of a rectangular Wilson loop of 
dimensions 2\r]\ x y/—£ 2 , see illustration Fig. 5.3b and the discussion in section 2.4.5 and 
Ref. [Col08]. The quantity 

B*™(P,Cft) = li{ ' £ ' vv> (5.26) 
VW(2M,V=F) 



4 Note that the third argument of the amplitudes also changes sign under combined (f) and (T). 
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does not need renormalization up to the quark field renormalization factor Z^. We plot the 
real part of i?^ vW in Fig. 5.4b. For short quark separations £, the divided ratio appears to 
reach a plateau value, just at the border to the gray shaded regions. According to eq. (5.25), 
this plateau value gives access to 2A 2 (divided by the Wilson loop factor): 

Re R^(P,C% V ) lr,v - Pl larse ; If w (£ 2 ,0,0,r\±l) , (5-27) 

where, for our parameters, £ 1//2 = \v-P\ = 0.52 GeV. Inside the shaded regions, the Wilson 
loop overlaps with itself on the periodic lattice (2|r/| > L = 20). Obviously, data points inside 
the shaded area must be discarded. For longer quark separations £, statistics breaks down 
before we can discern a plateau. However, we clearly observe that the divided amplitude 
^divW f lses w ith increasing quark separation. Whether this divided amplitude gives access to 
useful, well-defined TMD PDFs requires further studies. 



5.2.3 Ratios of Amplitudes, T-odd Effects from the Lattice 

Already in section 4.6.6 we noted that ratios of amplitudes are theoretically attractive, since 
they are inherently independent of the lattice renormalization scheme and scale. (That said, 
we should keep in mind that we found sizeable lattice cutoff effects for quark separations £ 
shorter than about 2a.) In the context of our test calculation with extended links, let us look 
at the quantity 

t> /p^iatw HP) Im R n( p A%) _ an + KWg) h 
° ddi ''' J - m N £- 2 P 1 ReR 7i (P,C^)- a 2 

±wP large A 12 (£ 2 , 0, 0, C\ ±1) + (m N I Pi) 2 B s (£ 2 , 0, 0, C 1 , ±1) (52g) 

I 2 (^ 2 ,o,o,c- 1 ,±i) 

where again C}/ 2 = 0.52 GeV, and where (m-Ar/Pj) 2 ~ 9.9 for the study at hand. The quantity 
above allows us to have a first glimpse at T-odd effects from the lattice. The signal we observe 
for i?odd contains a contribution from the Sivers function /^(^j Wl3b the amplitude A\ 2 , 
compare section 2.7.2. The values Ax 2 {- ■ ■ , 1) and • • , —1), differing only in sign, belong 
to different experimental setups, e.g., SIDIS versus Drell-Yan. 

Fig. 5.4c displays first encouraging results for i? dd- Firstly, the signal is clearly odd in rjvP, 
i.e., the lattice reproduces the expected T-odd effects. The data points for \/—£ 2 = 0.12 fm 
and 0.24 fm reach a stable plateau once \rjvP\ becomes greater than about 2. At a quark 
separation of 5a ~ 0.59 fm, it still looks like the data reaches a plateau value before statistics 
breaks down, and one may read off an estimate for the relative size of the linear combination 
A12 + {rriN /P\) 2 Bg compared to A\. 

Considering these results of our test calculation, it looks like it is possible to use lattice 
QCD for the calculation of "realistic" TMD PDFs, with gauge links as they occur in the 
description of experimental scattering processes such as SIDIS. Limitations are of technical, 
not of principal nature: the small nucleon momenta \ P\ available in present simulations restrict 
the evolution parameter £ to rather small values, and the proliferation of statistical noise as 
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we increase the staple extent rj confines us to small quark separations. Certainly, a smart 
choice of directions v, £, P and the spin projection matrices T 3pt is required to disentangle 
the 32 independent amplitudes Ai and Bi [GMS05] from the real and imaginary parts of the 
16 channels r op . As long as some open questions concerning the renormalization remain, we 
may resort to the calculation of ratios of amplitudes. Such ratios enable us to characterize 
the relative size of spin dependent phenomena (compare section 4.6.6) or, as in the present 
case, T-odd effects. 
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Figure 5.4: Results from a test calculation with staple-shaped link paths on the smeared 
coarse-m050 ensemble. Details regarding the setup are described in the text. 

(a) Real part of the unrenormalized ratio i? 7j (P, C\ a ^ v ). 

(b) Real part of the unrenormalized ratio divided by a Wilson loop subtraction factor. In the 
shaded region, the Wilson loop overlaps with its periodic image on the lattice. 

(c) Relative size of a linear combination of T-odd amplitudes to A2, obtained from a ratio of 
imaginary and real part of R Ti (P, C^L). 
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Conclusion 



6.1 Summary 

This work has been a first exploration of concepts, techniques, prospects and limitations 
for the calculation of transverse momentum dependent quark distribution functions (TMD 
PDFs) with lattice QCD. Comprehensive test calculations, mostly with a simplified operator 
(straight Wilson line), have produced encouraging first results and have enabled us to study 
renormalization properties of the relevant non-local operators in practice. We summarize the 
results as follows: 

• We can directly implement non-local operators q(£)U[Ce]q(0) on the lattice, by assem- 
bling the gauge link U[Cg] as a product of link variables connecting the two quark fields. 

• A gauge link of length I introduces a power divergence of the form exp(— 5m(a) I /a). 
Based on an analysis of smeared and unsmeared ensembles with four different lattice 
spacings, we confirm this behavior for our lattice operators with gauge links longer than 
about two lattice spacings. We have tested and compared methods to determine the 
renormalization constant dm. Two different non-perturbative methods, based on the 
static quark potential and on Wilson lines, prove to be successful and consistent with 
each other. The value Srh is defined unambiguously only with respect to a renormaliza- 
tion condition. 

• Fixing a renormalization condition, we are able to specify our results in terms of am- 
plitudes Ai unambiguously and independent of the lattice scheme and scale, up to a 
global multiplicative constant Z^\. The TMD PDFs are related to the corresponding 
amplitudes Ai by a Fourier transformation. 

• Restrictions from the Euclidean metric of the lattice preclude us from calculating the 
x-dependence of TMD PDFs directly. However, we do have access to Mellin moments. 
The first Mellin moment is just the x-integral of TMD PDFs and describes the k±- 
dependent distribution of quarks irrespective of their longitudinal momentum. 
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• From Gaussian fits to our amplitudes, we have calculated the first Mellin moment of 
selected TMD PDFs, namely fi 1)sW (k ± ), g^fix) and ^ sW (fc ± ). Note that: 

— We have connected the quark fields in our correlators with a straight Wilson line. 
The TMD PDFs we get (tagged with a superscript "sW" ) are therefore not strictly 
identical to those defined and used in the literature and for the description of, e.g., 
SIDIS. 

— In order to obtain unambiguous results, we need to fix a renormalization condition 
for 5rh. We choose a condition based on the static quark potential ((7 ren = 0, see 
section 4.4). The precise meaning of the renormalization condition in the context 
of factorization theorems for TMD PDFs has still to be worked out. 

— We work at a pion mass of about 500 MeV. 

— To be able to do the Fourier transform, we fit Gaussian functions to our amplitudes, 
excluding input data for quark separations smaller than 0.25 fm ~ 2a because of 
lattice cutoff effects. 

— The resulting Gaussian TMD PDFs clearly cannot reproduce the large-fcj_ behavior 
expected from perturbation theory, which is related to a singular dependence on £ 
at short distances. 

• We interpret our results in terms of quark densities as functions of the intrinsic trans- 
verse momentum of quarks inside the nucleon. The following observations are based on 
fi^ 8 ^(k±), ff[^ sW (fc±) an d dyr^ (k±) determined under the conditions and subject to 
the limitations mentioned above: 

— The unpolarized quark density in an unpolarized nucleon is axially symmetric in 
the transverse momentum plane. The width of the distribution, given by the root 
mean square transverse momentum, is (&j_)pi<xr = (381 ± 28) MeV for u — d quarks. 
This number is very sensitive to the renormalization condition we have imposed, 
but seems to be rather insensitive to the quark masses on our ensembles with pion 
masses ranging from about 500 to 800 MeV. 

— The density of longitudinally polarized quarks in a transversely polarized nu- 
cleon is no longer axially symmetric but appears deformed. This is due to a 
sizable contribution from <7^ sW . The peak of the density is shifted along the 
direction of the transverse nucleon spin vector. The direction of the shift is 
opposite for up and down quarks. The magnitude of the deformation can be 
characterized by the average transverse momentum, which we determine to be 
(k±)p TL = t( 73 ± 5 ) MeV l XS ± for U P quarks, and (k ± ) PTL = [(-31 ± 5) MeV] XS± 
for down quarks, all at a pion mass of about 500 MeV. These results are rather in- 
sensitive to the renormalization condition. Numbers compatible to the ones quoted 
above have been obtained with an alternative method, which involves extrapolation 
of a ratio of amplitudes that does not need renormalization. 

— As a cross check, we compute the axial vector coupling constant qa from the k±- 
integral of g^(k±), and find reasonable agreement with results in the literature 
extracted from the same ensembles. 
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• The dependence on the longitudinal momentum fraction x is encoded in the l-P- 
dependence of our amplitudes. Analyzing the renormalization scheme independent, nor- 
malized amplitude A 2 (£ 2 , £-P)/A 2 (£ 2 , 0), we observe factorization in £ 2 and l-P within 
our statistical and kinematical limits. This corresponds to the factorization /f w (x, k±) ~ 
fi(x)f[ 1)sW (k ± )/M. The results for our normalized amplitude are in qualitative agree- 
ment with a phenomenological parametrization of fi(x) and with a diquark model. 

• We carry out first studies with extended, staple-shaped gauge links, as they appear 
in the description of scattering experiments like SIDIS or the Drell-Yan process. The 
attainable evolution parameter £ = (P • v) 2 /\v 2 \ is limited by the maximum nucleon 
momentum \P\ achievable on the lattice. It is important to cancel the length dependent 
renormalization factors associated to the gauge link. For a time-reversal odd ratio of 
amplitudes, our test calculation yields a non-zero signal of encouraging quality. 

6.2 Open Questions and Future Projects 

Our study motivates a number of future research topics: 

• So far, we have studied only a small selection of TMD PDFs. A large number of other 
structures with interesting density interpretations can be readily analyzed with the 
techniques described in this work. 

• Can we detect the short distance behavior predicted by perturbation theory in our 
amplitudes obtained from the lattice? 

• In order to remove the power divergence associated to the gauge link, we have introduced 
a renormalization condition. How can we interpret this renormalization condition in 
terms of a continuum renormalization or factorization scale? 

• Can effective field theories or auxiliary field techniques give us hints regarding these 
issues? 

• It will be interesting to expand our test calculations with staple-shaped gauge links. 
Are the evolution equations in ( applicable to lattice results? Until all issues of renor- 
malization are resolved, one may resort to ratios of amplitudes, which should allow us 
to estimate the relative sizes of spin dependent or time-reversal-odd effects. 

• In particular in connection with "realistic", staple-shaped gauge links, how should we 
set up a Wilson loop subtraction factor to produce TMD PDFs that have a well-defined 
meaning as probability densities or as a non-perturbative ingredient to a factorized 
cross-section of a high-energy scattering process? 

• Can we improve our lattice operators by perturbatively motivated correction terms, 
similar as in the case of the static potential? 
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6.3 Resume 



Within this thesis, we have explored ways to calculate intrinsic transverse momentum distri- 
butions in hadrons with lattice QCD. There are remaining challenges to design well-defined 
observables that can describe these distributions. Using a simplified definition, we provide 
for the first time an insight into transverse momentum dependent quark distributions of the 
nucleon from a model independent calculation within full QCD. As an example, we deter- 
mine the strength of an axially asymmetric deformation of a spin-selective quark distribution 
inside the spin-polarized nucleon. Unlike most lattice calculations of hadron structure, we 
directly implement non-local operators and remove their length dependent divergence non- 
perturbatively. First test calculations indicate that lattice QCD has the potential to go 
beyond the simplified operator definition. It is well conceivable that the lattice perspective, 
being truly non-perturbative from the start and offering a numerical test bed, will be help- 
ful in developing a conceptually improved definition. Considering our initial success, we are 
optimistic that lattice QCD will become an important tool for the calculation of intrinsic 
transverse momentum dependent distributions of quarks. 



Appendix A 

Conventions and Useful Relations 



A.l Wilson Lines 

Let the continuous, piecewise differentiable function C(A), defined for < A < 1, parametrize 
a path from x to y, i.e., C XjJ ,(0) = x, C XtV (l) = y. We define a Wilson line (also called gauge 
link) along this path according to 

U\C\ ^Vex V (-igJ c MZ)) =Vexp(-ig J\\C'»(\) A„(C(A))) . (A.l) 

Here C'{\) = dC(\)/d\ is the derivative of the parametrization of the path. The path ordering 
symbol V indicates that the fields A fJi (C(X)) are sorted with increasing A from left to right 
after applying the definition of the matrix exponential. A Wilson line along a straight path 
from x to y will be denoted 

U[x, y\^V exp f-i g j" A^)\ . (A.2) 

Note that (U[x, y]y = U[y, x\. To keep notation concise, concatenations of Wilson lines along 
straight paths will be written as 

U[x,y]U[y,z]=U[x,y,z] (A.3) 
and accordingly for more than two line sections. 



A.2 Lightcone Coordinates 

We denote cartesian base vectors of Minkowski space no, • . ■ , h%- The non- vanishing compo- 
nents of the metric tensor g^ u = h^-hy are goo = — 9n = — 522 = —533 = 1- With base vectors 
"on the lightcone" given by 

n± = -7= (n ± n 3 ) , (A.4) 
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any four-vector v has a decomposition 

v = v + n + + w~n_ + v± , (A. 5) 

where v± is a spacelike vector in transverse direction, i.e. v± = v l h\ + v 2 fi2, v± • v± < 0. We 
also introduce a corresponding Euclidean two-component vector v± = (v 1 ,v 2 ) T , v± ■ v± > 0. 

In light cone notation, the metric tensor is given by = g |_ = —511 = —522 = 1- The 

scalar product of two general four-vectors v and w thus reads 

v • w = v + w ~ + v~w + + v± • w± = v + w~ + v~w + — v± • w± . (A. 6) 

Projections onto the transverse directions can be accomplished with the tensor g^ whose only 
non- vanishing components are —g±u = —g±22 = 1- The totally antisymmetric Levi-Civita 
symbol e tivpa follows the convention e 0123 = e 1-12 = 1. Its transverse projection e^" = e 
follows the convention e 12 = 1. In the two-dimensional, Euclidean transverse space we will 
denote it as a Levi-Civita symbol e±ij, with the transverse indices i,j € {1,2} and following 
the convention ej_ 12 = 1. 

We choose a nucleon momentum P on the light cone 

P = P+n + + . (A.7) 

In the frame we want to work in, the nucleon travels with large momentum in ^-direction, 
i.e. P + » m,N S> P~~ . The nucleon spin vector S shall be parameterized as 

S = -A ^ n_ + A — h + + S ± . (A.8) 
2P + mjv 

It fulfills S • P = and shall be normalized according to — S 2 = A 2 + S"j_ = 1. 



A. 3 Tensors in Minkowski Space 

e =+1, 7 = H 7 7 7 , ^ '7 J> CT 7 = ~jf <?a/3 ■ (A.9J 

A. 4 Dirac Spinors of free particles in Minkowski Space 

For a nucleon of momentum P and with a polarization vector 5 fulfilling S 2 = — 1, 5 • P = 0, 
we define the Dirac spinor J7 (P, S) such that 

17(P, S)U(P, S) = (f> + m N ) 1(1 - $^) . (A.10) 

Then 

1 = ^ U(P, S)U(P, S), S» = ^- U(P, S)j^U(P, S) . (A.ll) 



A. 5 Gordon Identities in Minkowski Space 



143 



A. 5 Gordon Identities in Minkowski Space 

_ _ 

U(P, S') 7 M U(P, S) = U(P, S') U(P, S) , (A. 12) 

m N 

U(P, S') 7^ 7 5 U(P, S) = U(P, S') ia^-f 5 — U(P, S) , (A. 13) 

m N 

U(P, S') a*"* U(P, S) = U(P, S') e^ a 7p7 5 — U(P, S) , (A. 14) 

m N 

U(P,S')^U(P,S)=0, (A.15) 
see, e.g., Ref. [DieOl]. 

A. 6 Trace Projections in Minkowski Space 

The relation between a Dirac matrix <]? and its projection $[ r ° p l is denned as 

$[r op ] = i trD ( $ r o P ) _ ( A l6 ) 

Then 

* = I 1 *M + I 7m $[7-1 + 1 ^ d>K1 _ 1 >7 5 ^] + 1 7 5 ^1 . (A . 17) 

A. 7 Euclidean Space 

Euclidean four- vectors shall be denoted xp, (p, = 1..4). We stick to common conventions: 

= ix^ = zxo j xj = = — Xj , ^1234 = "t - -^- • 

72 = 70 , 7J = -*y > 7S = 71727374 = ~7 5 , = |[7A> 7*] . # = 7£Pm • ( A - 18 ) 

where j = 1..3. It follows that 

d 3 = -ido , dj = dj . (A. 19) 

We have the following rules for rewriting an expression in Euclidean notation (in order of 
precedence): 

1. If we need to raise or lower an index before we can apply one of the rules below, we 
multiply by —1. 

2. ef^P^-ie^. 

3. - -a w , 7 5 - -75 • 
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4. 7^ -> % . 

5. d n . 

6. Upper indices \x become lower indices p,. 
Useful identities: 

{7p,,7v} = 18ji V , x-x = x^ = -XfiXfi = -x~x , 7^ = ij^dp , f=-ip. (A.20) 
The Dirac equation in Euclidean space reads 

($ + m)V> = 0. (A.21) 
The Dirac spinor of a free nucleon fulfills in Euclidean space 

U(P, S)U(P, S) = {-if + m N ) 1(1 - i$ Tb ) . (A.22) 
The trace projections keep their form: 

* = \l + \ 7, + \ * w ^ ] ~ \ 7,75 ^ ] + \ 75 <*> [7sl • (A.23) 
The Wilson line becomes 

U[C] ^Pexp(igJ c dtp = V exp (i g £ d\ C'-^X) A fi {C{\))\ . (A.24) 

For a path on the lattice C lat = (x^ n \x^ n ~ 1 \ . . . , x^ 2 \ x^\ 0) connecting the lattice sites at 
space-time locations x^\ we introduce gauge links as products of link variables according to 

w iat[ C iatj = U{x {n \x {n " l) ) ■■■ U(x (3 \x {2) )U(x {2 \x^)U(x^,0) . (A.25) 
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B.l Parametrization of the Correlators 



B.l.l M. structures 



M t =2A x , 

AV = — M P» + — A 3 
rriN rriN 



M 



A 4 {P fl k v -P v k^) 



l™ 2 N 



A 12 e^Pakp 7 „' 



+ A 9 e^Pp 7 a7 5 

m N 



+ — A w e^kp 7 «7 5 + \ A n ^k a Pp fa 5 , 



M 



m N 
-2 A 6 7^ 7 5 



m 



N 



M j5 = [2i A 5 h 5 ] ■ 



A 7 P^ 



A 8 k? H , 



Mi = 2 At 
2 



M 



m N 



A 2 P* x + 2im N A 3 P + 



2i A 12 e^PJp lv - 



M a n» = 



2% A 4 (P»t - P v e) + A 9 e^ a/3 P 7a7 5 

J m N 

.5 



+ 2im N Aio e^ a % 7 a7 5 -2m w A u e^ a H a P p / 7 5 , 
^ 7 m 7 s = -2 I 6 W - 2i A 7 P» / 7 5 + 2m N 2 A 8 t l~f> , 

-2m N A 5 / 7 5 1 



M 



An explicit expression for Mr°p is 

r d 4 £ 

M T o P (k,P;C e ) 



ik-i 



(2nY 



(0\^ N (P)q(£)r op U[C e } 9(0) A(P)|0) 
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Details 



where ip N (P) = f d 3 x e tP x ip N (x)\ x o =0 creates nucleons from the vacuum using the nucleon 
field iPn(x). 
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B.1.2 Table of Ratios 



QDP code 



T°p (Euclid.) 



r op (Mink.) 



o 
l 

2 
3 
4 
5 
6 
7 
8 
9 
10 



11 
12 
13 
14 
15 



1 

71 
75 
2 bi,72] 
73 

I [71 ,75] 

I [72 ,73] 
-727S 
72 

I [71, 72] 

1 [72 >72] 

7375 

2 [73 ,72] 
-727S 

7175 
75 



-i7 X 

-i7 2 
12 



IG 



1(7 



Z7 

13 



23 



10 



7 7 5 



7 



01 



02 



a 



03 



a 



-i7 2 7 5 



■ 1 ^ 

vy 7 



-7 



P(P) 

i 



P(P) 



A 2 Pi + 



JV 



P(P) 



m 



N 



E{P) 
m N 
E{P) 



A 3 £ 2 + tun Ai 2 £i 



I 4 ^ 2 Pi + i A 9 + im 2 N i n (£3? 



in 



N 



E{P) 
m N 



A3 £3 



E(P) 



A i £ 3 P 1 -im 2 N A u £2£3 



im N A u £±£3 
irriN A7 £3 
A 2 



mN A l2 £ 2Pl 



E{P) 



im 



im N A A £1 - A w £ 2 

E{P) 



irriN A^ £ 2 + 



1 



E{P) 



A gPl 



im 



v A 10 £ 1 + ^An(£3) 2 Pi 



E(P) 



E{P) 



A. 



im 



N 



E(P) 
ivnjq A4 £3 



E(P) 



A 



8 1*3 J 



N 



E(P) 



A ll £ 2 £ z P 1 



vm"* 



E(P) 823 



im 



N 



E{P) 



m: 



N 



A s £i£ 3 



A 5 £3 



rnN 
E(P) 



A 7 £ 3 Pi 



Table B.l: Plateau values of the ratios Pp e p(T, P,£) in terms of the amplitudes A4. Here we 
employ the LHPC conventions for T 2pt and r 3pt , i.e. the nucleons are spin-projected along 
the z-axis. We choose the nucleon momentum P = (Pi, 0,0) = (Pi, 0,0), and the quark 
separation is £ = (£i,£ 2 ,£ 3 ) = (£i,£ 2 ,£ 3 ). Naturally, we have £° = £\ = on the lattice. We 
scan through the different Dirac contractions r op , sticking to the conventions of the QDP++ 
programming library [Edw07]. Note that A4, A§ and Ai 2 vanish for straight link paths. 
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B.l. 3 General Fourier Transformation relating the Ai to TMD PDFs 

From eq. (4.14), we see that we need to rewrite Fourier transformations of the form 



Y 



31 ,32, — ,3n 



d(e-p) i(e . P)x fd 2 £ ± Ux . kj 

27T J (27T)2 



x£ ±j e ±j2 ---e ±jn (t) m F(f,e-p) 



l+=0 



as an integration with respect to Lorentz-invariant quantities. Note that £ 2 
first concentrate on the integral with respect to £±: 

d 2 e 



e+=o 



Y 



(2vr) s 



i£± ■ k± 



£±n i± j2 ---£± jn F(-£i,£-P) 



(B.4) 
We 

(B.5) 



We express the i±j in terms of derivatives and integrate with respect to the angular degree 
of freedom: 



Y± ■ ■ 

31,32, -,3n 







dk ±jl 







d 2 l. 



At i k 



d 



dk ±j J J (2vr)2 

"°° d(-£ 2 ) f 2w do 



F(-£ 2 ± ,£-P) 



dk ±jl 



) { l dk ±j J J 2(2vr) 7o 2vr H ^ F) 











d(-e 2 ) 



Jo(y/^\k ± \) F(£ 2 ,£-P) . (B.6) 



dk ±jl J \ dk ±jn J J 2(2tt) 
Obviously the integral only depends on the absolute value of k±, so we can replace 
d d\k±\ d kj_j i — — d 



dk±j dk±j d\k±\ \kj 



d[V^P\k_ 



rriN z oz rriN 



(B.7) 



In the last step we have introduced the abbreviation z = y—£? | k± \ . We now make use of a 
property of the Bessel function [RWVOO] 



I d ( J p {z)\ J p +i(z) 



zdz \ zP J zP +1 
which holds for any non-negative integer p. Applying this iteratively in (B.6), we get 



Y 1 - ■ ■ = i 

31,32, -,3n 



n K ±3l K ±3n 



N 



d(-£ 2 ) J n (V^\k ± \) 

o 2 ( 27 (V=F|fe ± |) n 



(-t 2 m 2 N r 



1 



m 



F(£ 2 ,£-P) 



N 



(B.9) 

Regarding the integral with respect to I ■ P, we only have to substitute I = (£ ■ P)/P + \g +=Q 
in eq. (B.4). Summarizing our results, we can write 



Y 



k±j 1 ■ ■ ■ k±j n 



3l,32,-,3n 



in 



T m ,n{x, \k_\ 



N 



1 

(P+) m m 



F(£ 2 ,£-P) 



N 



(B.10) 
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where we have defined 



T m , n (x, \k±\) [(expression)] = 



.„ f d(£_P) e _ t(e . P)x {e _ p)m /°° d ^f} Y^ lk± , l } (-fm 2 N r (expression) . (B.ll) 



o 2(2vr) (^\ k 



2vr 



Note that the operation T m ^ n has the dimension mass 2 . 



B.2 Ingredients to the Perturbative Calculation 

The inverse gluon propagator for the MILC gluon action is given by [Bis05] 

Dfiv(k) = k^kp + g^p ^(l - (c 2 + c 3 )k 2 ) y]«f - k 4 p - kpkpk 2 ^ (cl - c2 - c3) ^ 



p(l - gpp) ( 1 - Cl(«J + K?) - (C2 + c 3 ) £ /c 2 , ] , (B.12) 



where kn = 2 sin(akn/2). Note that the expression Dnp(k) alone corresponds to the gauge 
with £ = 1 in eq. (4.41). The parameters c, are determined according to 

c = 5/3, ci = — ^(1 + 0.4805a s ) , c 2 = , c 3 = ^0.03325a s . (B.13) 
12uq 3uq 

where a s = — 4 log (no )/3. 0684. 

For perturbative calculations with HYP blocking, the coefficients hpp(k) are given by [DeG03, 
Bis05] 

hp,p{k) = 5p_p ^1 — ^ ^ kp J7^p(A:)^ H q^^P Ki> ^P&(.ty ' 

n A p(fc) = 1 + q 2 (1 + a 3 ) - ^(1 + 2a 3 )(k 2 - k\ - kl) + J] k % ■ (B-14) 

where the coefficients a« have been chosen as specified in section 3.1.15. 



B.3 Gluon-Exchange Corrections to the Static Potential 

In Fig. B.l we show the size and shape of the corrections to the static quark potential as we 
apply them in our fits in section 4.4.7.2. We remark that the large uncertainty in A for the 
unsmeared ensemble is a harmless symptom caused by the small number of configurations 
and a reduced set of lattice vectors r evaluated for the unsmeared ensemble. The qualitative 
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Figure B.l: Corrections to the static quark potential as obtained from lattice perturbation 
theory for the coarse- m020 ensemble, plotted versus f = \r\ for various lattice vectors r . For 
A, we have taken the central value from our fits to lattice data. 

(a) Unsmeared ensemble, A = 0.30(35), 

(b) Smeared ensemble, A = 0.442(29). 



features of the corrections plotted in Fig. B.la resemble very much those shown for a different 
action in Ref. [BalOl], in particular we note that the sign of the corrections alternates quickly. 
In contrast, for the HYP smeared ensemble, the corrections are systematically negative, and 
much larger in size. Provided our implementation of the HYP smeared propagator is correct, 
these results are an indication that HYP smearing introduces sizeable lattice artefacts at 
distances smaller than about three lattice spacings. The displeasing feature of these artefacts 
is that they are not easy to recognize in the data because of their "smooth" dependence on r. 



B.4 Fourier Transform of the High-fc^ Behavior 

Suppose ^ + ^\k ± ;P,S) = fP(k±) « b/k± for some real constant b above some threshold 
\k±\ > &min- The corresponding contribution to the amplitude A2(£ 2 ,0) reads 

J\k ± \>k min k ± Jk min \K±\ 

= 2nb[" du J ^ + 2, b r du J ^. (B.15) 

We now make an approximation valid for A; m i n |£j_| <k<1. Taylor expansion yields Jo(u) = 
1 + 0(u 2 ), so we get 

27Tb In f — -j— t\ + 2nb f°° du + (k 2 ) + 0(k 2 min \e ± \ 2 ) = 

- 27Tb In (M\ + O {(k min \£ ± \) 2 ) . (B.16) 

for an appropriate choice of Iq. 
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B.5 Error Propagation for the Gaussian Parametrization 

Let us derive a simple estimate how systematic errors in 5m affect the width of our Gaussian 
fits. Suppose we "fit" a Gaussian curve to the lattice data at just two quark separations l\ 
and I2: 



c i expHf/a]) = 2i- 1 2 exp(-6mh) i™(-Z?,0) , 
c,-expH!/a]) = ^ exp(-5mZ 2 ) i™(-Z 2 2 ,0) . 



(B.17) 



Solving this for l/c| and differentiating with respect to 5m, we obtain a formula for the error 
in the inverse width of the Gaussian: 







d(5m) 



h + h 



A 



07 



a(h + Z 2 



-A[5m] 



(B.18) 



For l\ and I2, we simply take our fit range: l\ = 0.25 fm, / 2 = -L/2. 



B.6 k -Moments 

In general, fc^-moments of the TMD quark-quark correlator are of the form 

J d 2 k ± k ±n ■ ■ ■ k ±jn ^ r ° P ^\k ± ; P, S) . (B.19) 

Assuming the integral above exists, and assuming that ^ r ° P \i, P, S) vanishes quickly enough 
for large £±, we may we may rewrite the expression above with the help of eq. (2.13) as 



After n partial integrations, and carrying out the Fourier transform, we obtain 

^a£--a£-* ,r "'(^' s »l-- < B ' 21 > 

In order to analyze a particular channel, we insert the parametrization of $[ r ° p l in terms of 
amplitudes Ai according to eq. (4.5). Here we can make use of the relation 

^- Mf,0) = 2£ ±j ^—^(^,0) . (B.22) 



B.6 fcj_-Moments 
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Thus we find 



J d 2 k ± ^ 1 \k ± ;P,S) = 2A 2 (0,0) , 
J d 2 k ± k ± ^ 1 \k r ,P,S) = 2i£ ± ^^2A 2 (i 2 ,0)\ e=0 , 

J d 2 k ± kl ^ +](1) (^;P,5) = -4(^ + M 2 )^^) 2 I 2 (£ 2 ,o)| £=0 

J d 2 k ± $ [7+75l(1) (fc ± ;P,5)| s±=0 = -2Al 6 (0,0) , 
J d 2 k ± k ± ^ 5]W (k ± ;P,S)\ A=Q = - (s ± + 2£ ± (S ± -£ ± )^^j 2m N A 7 (£ 2 ,0) 



i=o ■ 
(B.23) 
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